!> \file cu_ntiedtke.F90
!! This file contains the CCPP-compliant new Tiedtke scheme which parameterize
!! Shallow, deep, and mid-level convections in the model
!! Please refer to Tiedtke (1989), Bechtold et al. (2004,2008, 2014),
!! Zhang et al. (2011), Zhang and Wang (2017, 2018)
!!
!###########################################################

module cu_ntiedtke

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     use machine , only : kind_phys
     ! DH* TODO - replace with arguments to subroutine calls,
     ! this also requires redefining derived constants in the
     ! parameter section below
     use physcons, only:rd=>con_rd, rv=>con_rv, g=>con_g, &
   &       cpd=>con_cp, alv=>con_hvap, alf=>con_hfus              

     implicit none
     real(kind=kind_phys),private :: rcpd,vtmpc1,tmelt,als,t13,                &
             c1es,c2es,c3les,c3ies,c4les,c4ies,c5les,c5ies,zrg

     real(kind=kind_phys),private :: rovcp,r5alvcp,r5alscp,ralvdcp,ralsdcp,ralfdcp,rtwat,rtber,rtice 
     real(kind=kind_phys),private :: entrdd,cmfcmax,cmfcmin,cmfdeps,zdnoprc,cprcon
     integer,private :: momtrans,p650

     parameter(         &
      t13 = 0.333333333,&
      rcpd=1.0/cpd,     &
      tmelt=273.16,     &
      zrg=1.0/g,        &
      c1es=610.78,      &
      c2es=c1es*rd/rv,  &
      c3les=17.2693882, &
      c3ies=21.875,     &
      c4les=35.86,      &
      c4ies=7.66,       &
      als  = alv+alf,   &
      c5les=c3les*(tmelt-c4les),  &
      c5ies=c3ies*(tmelt-c4ies),  &
      r5alvcp=c5les*alv*rcpd,     &
      r5alscp=c5ies*als*rcpd,     &
      ralvdcp=alv*rcpd,           &
      ralsdcp=als*rcpd,           &
      ralfdcp=alf*rcpd,           &
      rtwat=tmelt,                &
      rtber=tmelt-5.,             &
      rtice=tmelt-23.,            &
      vtmpc1=rv/rd-1.0,           &
      rovcp = rd*rcpd )
!
!     entrdd: average entrainment & detrainment rate for downdrafts
!     ------
!
      parameter(entrdd = 2.0e-4)
!
!     cmfcmax:   maximum massflux value allowed for updrafts etc
!     -------
!
      parameter(cmfcmax = 1.0)
!
!     cmfcmin:   minimum massflux value (for safety)
!     -------
!
      parameter(cmfcmin = 1.e-10)
!
!     cmfdeps:   fractional massflux for downdrafts at lfs
!     -------
!
      parameter(cmfdeps = 0.30)

!     zdnoprc:   deep cloud is thicker than this height (Unit:Pa)
!
      parameter(zdnoprc = 2.0e4)
!     -------
!   
!     cprcon:    coefficient from cloud water to rain water
!   
      parameter(cprcon = 1.4e-3)
!     -------
!
!     momtrans: momentum transport method
!     ( 1 = IFS40r1 method; 2 = new method )
!
      parameter(momtrans = 2 )
!     -------
!
      logical :: isequil
!     isequil: representing equilibrium and nonequilibrium convection
!     ( .false. [default]; .true. [experimental]. Ref. Bechtold et al. 2014 JAS )
! 
      parameter(isequil = .false. )
!
!--------------------
!     switches for deep, mid, shallow convections, downdraft, and momemtum transport
!     ------------------
      logical :: lmfpen,lmfmid,lmfscv,lmfdd,lmfdudv
      parameter(lmfpen=.true.,lmfmid=.true.,lmfscv=.true.,lmfdd=.true.,lmfdudv=.true.)
!--------------------
!#################### end of variables definition##########################
!-----------------------------------------------------------------------
!
contains
!> \brief Brief description of the subroutine
!!
!! \section arg_table_cu_ntiedtke_init Argument Table
!! \htmlinclude cu_ntiedtke_init.html
!!
      subroutine cu_ntiedtke_init(imfshalcnv, imfshalcnv_ntiedtke, imfdeepcnv,  &
                          imfdeepcnv_ntiedtke,mpirank, mpiroot, errmsg, errflg)

         implicit none

         integer,                   intent(in) :: imfshalcnv, imfshalcnv_ntiedtke
         integer,                   intent(in) :: imfdeepcnv, imfdeepcnv_ntiedtke           
         integer,                   intent(in)    :: mpirank
         integer,                   intent(in)    :: mpiroot
         character(len=*),          intent(  out) :: errmsg
         integer,                   intent(  out) :: errflg
         
         ! initialize ccpp error handling variables
         errmsg = ''
         errflg = 0
         
         ! DH* temporary
         if (mpirank==mpiroot) then
            write(0,*) ' -----------------------------------------------------------------------------------------------------------------------------'
            write(0,*) ' --- WARNING --- the CCPP New Tiedtke convection scheme is currently under development, use at your own risk --- WARNING ---'
            write(0,*) ' -----------------------------------------------------------------------------------------------------------------------------'
         end if
         ! *DH temporary

         ! Consistency checks
         if (imfshalcnv/=imfshalcnv_ntiedtke) then
           write(errmsg,'(*(a))') 'Logic error: namelist choice of',       &
        &    ' shallow convection is different from new Tiedtke scheme'
           errflg = 1
           return
         end if

         if (imfdeepcnv/=imfdeepcnv_ntiedtke) then
           write(errmsg,'(*(a))') 'Logic error: namelist choice of',       &
        &    ' deep convection is different from new Tiedtke scheme'
           errflg = 1
           return
         end if
         
      end subroutine cu_ntiedtke_init

! Tiedtke cumulus scheme from WRF with small modifications
! This scheme includes both deep and shallow convections
!===================
!
!> \section arg_table_cu_ntiedtke_run Argument Table
!! \htmlinclude cu_ntiedtke_run.html
!!
!-----------------------------------------------------------------------
!          level 1 subroutine 'tiecnvn'
!-----------------------------------------------------------------
      subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, &
                                 evap,hfx,zprecc,lmask,lq,km,dt,dx,kbot,ktop,kcnv, &
                                 ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,errmsg,errflg)
!-----------------------------------------------------------------
!  this is the interface between the model and the mass 
!  flux convection module
!-----------------------------------------------------------------
      implicit none
! in&out variables
      integer, intent(in)  ::  lq, km, ktrac
      real(kind=kind_phys),     intent(in ) :: dt
      integer, dimension( : ),   intent(in)  :: lmask
      real(kind=kind_phys), dimension( : ),     intent(in ) :: evap, hfx, dx
      real(kind=kind_phys), dimension( :, : ),     intent(inout) :: pu, pv, pt, pqv
      real(kind=kind_phys), dimension( :, :),     intent(in )   :: tdi, qvdi, poz, prsl, pomg, pqvf, ptf
      real(kind=kind_phys), dimension( :, : ),   intent(in )   :: pzz, prsi
      real(kind=kind_phys), dimension( :, :, : ),    intent(inout ) ::  clw

      integer, dimension( : ),   intent(out)  :: kbot, ktop, kcnv
      real(kind=kind_phys), dimension( : ),   intent(out)  :: zprecc
      real(kind=kind_phys), dimension (:, :), intent(out)  :: ud_mf, dd_mf, dt_mf, cnvw, cnvc
 
! error messages
      character(len=*), intent(out)    ::                                 errmsg
      integer,          intent(out)    ::                                 errflg
 
! local variables
      real(kind=kind_phys) pum1(lq,km),  pvm1(lq,km),  ztt(lq,km),                &
     &     ptte(lq,km),    pqte(lq,km),  pvom(lq,km),  pvol(lq,km),               &
     &     pverv(lq,km),   pgeo(lq,km),  pap(lq,km),   paph(lq,km+1)
      real(kind=kind_phys) pqhfl(lq),      zqq(lq,km),                            &
     &     prsfc(lq),      pssfc(lq),    pcte(lq,km),                             &
     &     phhfl(lq),      pgeoh(lq,km+1)
      real(kind=kind_phys) ztp1(lq,km),    zqp1(lq,km),  ztu(lq,km),   zqu(lq,km),&
     &     zlu(lq,km),     zlude(lq,km), zmfu(lq,km),  zmfd(lq,km),  zmfude_rate(lq,km),&
     &     zqsat(lq,km),   zrain(lq)
      real(kind=kind_phys),allocatable ::  pcen(:,:,:),ptenc(:,:,:) 

      integer icbot(lq),   ictop(lq),     ktype(lq),   lndj(lq)
      logical locum(lq)
!
      real(kind=kind_phys) ztmst,fliq,fice,ztc,zalf,tt
      integer i,j,k,k1,n,km1,ktracer
      real(kind=kind_phys) ztpp1
      real(kind=kind_phys) zew,zqs,zcor
!
! Initialize CCPP error handling variables
   errmsg = ''
   errflg = 0

      km1 = km + 1
      ztmst=dt
!
!  masv flux diagnostics.
!
      do j=1,lq
        zrain(j)=0.0
        locum(j)=.false.
        prsfc(j)=0.0
        pssfc(j)=0.0
        pqhfl(j)=evap(j)
        phhfl(j)=hfx(j)
        pgeoh(j,km1)=pzz(j,1)
         paph(j,km1)=prsi(j,1)
        if(lmask(j).eq.1) then
           lndj(j)=1
        else
           lndj(j)=0
        end if
      end do
!
!     convert model variables for mflux scheme
!
      do k=1,km
        k1=km-k+1
        do j=1,lq
          pcte(j,k1)=0.0
          pvom(j,k1)=0.0
          pvol(j,k1)=0.0
          ztp1(j,k1)=pt(j,k)
          zqp1(j,k1)=pqv(j,k)
          pum1(j,k1)=pu(j,k)
          pvm1(j,k1)=pv(j,k)
          pverv(j,k1)=pomg(j,k)
          pgeo(j,k1)=poz(j,k)
          pgeoh(j,k1)=pzz(j,k+1)
          pap(j,k1)=prsl(j,k)
          paph(j,k1)=prsi(j,k+1)
          tt=ztp1(j,k1)
          zew  = foeewm(tt)                                      
          zqs  = zew/pap(j,k1)                                      
          zqs  = min(0.5,zqs)                                         
          zcor = 1./(1.-vtmpc1*zqs)                                    
          zqsat(j,k1)=zqs*zcor      
          pqte(j,k1)=pqvf(j,k)+(pqv(j,k)-qvdi(j,k))/ztmst 
          zqq(j,k1) =pqte(j,k1)
          ptte(j,k1)=ptf(j,k)+(pt(j,k)-tdi(j,k))/ztmst
          ztt(j,k1) =ptte(j,k1)
          ud_mf(j,k1)=0.
          dd_mf(j,k1)=0.
          dt_mf(j,k1)=0.
          cnvw(j,k1)=0.
          cnvc(j,k1)=0.
        end do
      end do

      if(ktrac > 2) then
        ktracer = ktrac - 2
        allocate(pcen(lq,km,ktracer))
        allocate(ptenc(lq,km,ktracer))
        do n=1,ktracer
          do k=1,km
            k1=km-k+1
            do j=1,lq
              pcen(j,k1,n) = clw(j,k,n+2)
              ptenc(j,k1,n)= 0.
            end do
          end do
        end do
      else
        ktracer = 2
        allocate(pcen(lq,km,ktracer))
        allocate(ptenc(lq,km,ktracer))
        do n=1,ktracer
          do k=1,km
            do j=1,lq
              pcen(j,k,n) = 0.
              ptenc(j,k,n)= 0.
            end do
          end do
        end do
      end if
      
!      print *, "pgeo=",pgeo(1,:)
!      print *, "pgeoh=",pgeoh(1,:)
!      print *, "pap=",pap(1,:)
!      print *, "paph=",paph(1,:)
!      print *, "ztp1=",ztp1(1,:)
!      print *, "zqp1=",zqp1(1,:) 
!      print *, "pum1=",pum1(1,:)
!      print *, "pvm1=",pvm1(1,:)
!      print *, "pverv=",pverv(1,:)
!      print *, "pqte=",pqte(1,:)
!      print *, "ptte=",ptte(1,:)
!      print *, "hfx=", pqhfl(1),phhfl(1),dx(1)
!
!-----------------------------------------------------------------------
!*    2.     call 'cumastrn'(master-routine for cumulus parameterization)
!
      call cumastrn        &
     &    (lq,       km,       km1,      km-1,    ztp1, &
     &     zqp1,     pum1,     pvm1,     pverv,   zqsat,&
     &     pqhfl,    ztmst,    pap,      paph,    pgeo, &
     &     ptte,     pqte,     pvom,     pvol,    prsfc,&
     &     pssfc,    locum,    ktracer,  pcen,    ptenc,&
     &     ktype,    icbot,    ictop,    ztu,     zqu,  &
     &     zlu,      zlude,    zmfu,     zmfd,    zrain,&
     &     pcte,     phhfl,    lndj,     pgeoh,   zmfude_rate, dx)
!
!     to include the cloud water and cloud ice detrained from convection
!
      do k=1,km
      k1=km-k+1
      do j=1,lq
      if(pcte(j,k1).gt.0.) then
        fliq=foealfa(ztp1(j,k1))
        fice=1.0-fliq
        clw(j,k,2)=clw(j,k,2)+fliq*pcte(j,k1)*ztmst
        clw(j,k,1)=clw(j,k,1)+fice*pcte(j,k1)*ztmst
      endif
      end do
      end do
!
      do k=1,km
        k1 = km-k+1
        do j=1,lq
          pt(j,k) = ztp1(j,k1)+(ptte(j,k1)-ztt(j,k1))*ztmst
          pqv(j,k)= zqp1(j,k1)+(pqte(j,k1)-zqq(j,k1))*ztmst
          ud_mf(j,k)= zmfu(j,k1)*ztmst
          dd_mf(j,k)= -zmfd(j,k1)*ztmst
          dt_mf(j,k)= zmfude_rate(j,k1)*ztmst
          cnvw(j,k) = zlude(j,k1)*ztmst*g/(prsi(j,k)-prsi(j,k+1))
          cnvc(j,k) = 0.04 * log(1. + 675. * ud_mf(j,k))
          cnvc(j,k) = min(cnvc(j,k), 0.6)
          cnvc(j,k) = max(cnvc(j,k), 0.0)
        end do
      end do

      do j=1,lq
        zprecc(j)=amax1(0.0,(prsfc(j)+pssfc(j))*ztmst*0.001)
        kbot(j) = km-icbot(j)+1
        ktop(j) = km-ictop(j)+1
        if(ktype(j).eq.1 .or. ktype(j).eq.3) then
           kcnv(j)=1 
        else
           kcnv(j)=0
        end if
      end do

      if (lmfdudv) then
        do k=1,km
          k1=km-k+1
          do j=1,lq
            pu(j,k)=pu(j,k)+pvom(j,k1)*ztmst
            pv(j,k)=pv(j,k)+pvol(j,k1)*ztmst
          end do
        end do
      endif

!
! Currently, vertical mixing of tracers are turned off
!      if(ktrac > 2) then
!        do n=1,ktrac-2
!          do k=1,km
!            k1=km-k+1
!            do j=1,lq
!              clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst
!            end do
!          end do
!        end do
!      end if
      deallocate(pcen)
      deallocate(ptenc)
!
      return
      end subroutine cu_ntiedtke_run

!#############################################################
!
!             level 2 subroutines
!
!#############################################################
!***********************************************************
!           subroutine cumastrn
!***********************************************************
      subroutine cumastrn  &
     &    (klon,     klev,     klevp1,   klevm1,   pten, &
     &     pqen,     puen,     pven,     pverv,    pqsen,&
     &     pqhfl,    ztmst,    pap,      paph,     pgeo, &
     &     ptte,     pqte,     pvom,     pvol,     prsfc,& 
     &     pssfc,    ldcum,    ktrac,    pcen,     ptenc,&
     &     ktype,    kcbot,    kctop,    ptu,      pqu,&
     &     plu,      plude,    pmfu,     pmfd,     prain,&
     &     pcte,     phhfl,    lndj,     zgeoh,    pmfude_rate, dx)
      implicit none
!
!***cumastrn*  master routine for cumulus massflux-scheme
!     m.tiedtke      e.c.m.w.f.      1986/1987/1989
!     modifications
!     y.wang         i.p.r.c         2001
!     c.zhang                        2012
!***purpose
!   -------
!          this routine computes the physical tendencies of the
!     prognostic variables t,q,u and v due to convective processes.
!     processes considered are: convective fluxes, formation of
!     precipitation, evaporation of falling rain below cloud base,
!     saturated cumulus downdrafts.
!***method
!   ------
!     parameterization is done using a massflux-scheme.
!        (1) define constants and parameters
!        (2) specify values (t,q,qs...) at half levels and
!            initialize updraft- and downdraft-values in 'cuinin'
!        (3) calculate cloud base in 'cutypen', calculate cloud types in cutypen,
!            and specify cloud base massflux
!        (4) do cloud ascent in 'cuascn' in absence of downdrafts
!        (5) do downdraft calculations:
!              (a) determine values at lfs in 'cudlfsn'
!              (b) determine moist descent in 'cuddrafn'
!              (c) recalculate cloud base massflux considering the
!                  effect of cu-downdrafts
!        (6) do final adjusments to convective fluxes in 'cuflxn',
!            do evaporation in subcloud layer
!        (7) calculate increments of t and q in 'cudtdqn'
!        (8) calculate increments of u and v in 'cududvn'
!***externals.
!   ----------
!       cuinin:  initializes values at vertical grid used in cu-parametr.
!       cutypen: cloud bypes, 1: deep cumulus 2: shallow cumulus
!       cuascn:  cloud ascent for entraining plume
!       cudlfsn: determines values at lfs for downdrafts
!       cuddrafn:does moist descent for cumulus downdrafts
!       cuflxn:  final adjustments to convective fluxes (also in pbl)
!       cudqdtn: updates tendencies for t and q
!       cududvn: updates tendencies for u and v
!***switches.
!   --------
!          lmfmid=.t.   midlevel convection is switched on
!          lmfdd=.t.    cumulus downdrafts switched on
!          lmfdudv=.t.  cumulus friction switched on
!***
!     model parameters (defined in subroutine cuparam)
!     ------------------------------------------------
!     entrdd     entrainment rate for cumulus downdrafts
!     cmfcmax    maximum massflux value allowed for
!     cmfcmin    minimum massflux value (for safety)
!     cmfdeps    fractional massflux for downdrafts at lfs
!     cprcon     coefficient for conversion from cloud water to rain
!***reference.
!   ----------
!          paper on massflux scheme (tiedtke,1989)
!-----------------------------------------------------------------
      integer  klev,klon,ktrac,klevp1,klevm1
      real(kind=kind_phys)     pten(klon,klev),        pqen(klon,klev),& 
     &         puen(klon,klev),        pven(klon,klev),&
     &         ptte(klon,klev),        pqte(klon,klev),&
     &         pvom(klon,klev),        pvol(klon,klev),&
     &         pqsen(klon,klev),       pgeo(klon,klev),&
     &         pap(klon,klev),         paph(klon,klevp1),&
     &         pverv(klon,klev),       pqhfl(klon),&
     &         phhfl(klon)            
      real(kind=kind_phys)     ptu(klon,klev),         pqu(klon,klev),&
     &         plu(klon,klev),         plude(klon,klev),&
     &         pmfu(klon,klev),        pmfd(klon,klev),&
     &         prain(klon),&
     &         prsfc(klon),            pssfc(klon)
      real(kind=kind_phys)     ztenh(klon,klev),       zqenh(klon,klev),&
     &         zgeoh(klon,klevp1),     zqsenh(klon,klev),&
     &         ztd(klon,klev),         zqd(klon,klev),&
     &         zmfus(klon,klev),       zmfds(klon,klev),&
     &         zmfuq(klon,klev),       zmfdq(klon,klev),&
     &         zdmfup(klon,klev),      zdmfdp(klon,klev),&
     &         zmful(klon,klev),       zrfl(klon),&
     &         zuu(klon,klev),         zvu(klon,klev),&
     &         zud(klon,klev),         zvd(klon,klev),&
     &         zlglac(klon,klev)
      real(kind=kind_phys)     pmflxr(klon,klevp1),    pmflxs(klon,klevp1)
      real(kind=kind_phys)     zhcbase(klon),& 
     &         zmfub(klon),            zmfub1(klon),&
     &         zdhpbl(klon)          
      real(kind=kind_phys)     zsfl(klon),             zdpmel(klon,klev),&
     &         pcte(klon,klev),        zcape(klon),&
     &         zcape1(klon),           zcape2(klon),&
     &         ztauc(klon),            ztaubl(klon),&
     &         zheat(klon)
      real(kind=kind_phys)     pcen(klon,klev,ktrac), ptenc(klon,klev,ktrac) 
      real(kind=kind_phys)     wup(klon),              zdqcv(klon)            
      real(kind=kind_phys)     wbase(klon),            zmfuub(klon)
      real(kind=kind_phys)     upbl(klon)
      real(kind=kind_phys)     dx(klon)
      real(kind=kind_phys)     pmfude_rate(klon,klev), pmfdde_rate(klon,klev)
      real(kind=kind_phys)     zmfuus(klon,klev),      zmfdus(klon,klev)
      real(kind=kind_phys)     zmfudr(klon,klev),      zmfddr(klon,klev) 
      real(kind=kind_phys)     zuv2(klon,klev),ztenu(klon,klev),ztenv(klon,klev)
      real(kind=kind_phys)     zmfuvb(klon),zsum12(klon),zsum22(klon)
      integer  ilab(klon,klev),        idtop(klon),&
     &         ictop0(klon),           ilwmin(klon)
      integer  kdpl(klon)
      integer  kcbot(klon),            kctop(klon),&
     &         ktype(klon),            lndj(klon)
      logical  ldcum(klon),            lldcum(klon)
      logical  loddraf(klon),          llddraf3(klon), llo1,   llo2(klon)

!  local varaiables
      real(kind=kind_phys)     zcons,zcons2,zqumqe,zdqmin,zdh,zmfmax
      real(kind=kind_phys)     zalfaw,zalv,zqalv,zc5ldcp,zc4les,zhsat,zgam,zzz,zhhat
      real(kind=kind_phys)     zpbmpt,zro,zdz,zdp,zeps,zfac,wspeed
      integer  jl,jk,ik
      integer  ikb,ikt,icum,itopm2
      real(kind=kind_phys)     ztmst,ztau,zerate,zderate,zmfa
      real(kind=kind_phys)     zmfs(klon),pmean(klev),zlon
      real(kind=kind_phys)     zduten,zdvten,ztdis,pgf_u,pgf_v
!-------------------------------------------
!     1.    specify constants and parameters
!-------------------------------------------
      zcons=1./(g*ztmst)
      zcons2=3./(g*ztmst)

      zlon = real(klon)
      do jk = klev , 1 , -1
        pmean(jk) = sum(pap(:,jk))/zlon
      end do
      p650 = klev-2
      do jk = klev , 3 , -1
        if ( pmean(jk)/pmean(klev)*1.013250e5 >  650.e2 ) p650 = jk
      end do

!--------------------------------------------------------------
!*    2.    initialize values at vertical grid points in 'cuini'
!--------------------------------------------------------------
      call cuinin &
     &    (klon,     klev,     klevp1,   klevm1,   pten, &
     &     pqen,     pqsen,    puen,     pven,     pverv,&
     &     pgeo,     paph,     zgeoh,    ztenh,    zqenh,&
     &     zqsenh,   ilwmin,   ptu,      pqu,      ztd,  &
     &     zqd,      zuu,      zvu,      zud,      zvd,  &
     &     pmfu,     pmfd,     zmfus,    zmfds,    zmfuq,&
     &     zmfdq,    zdmfup,   zdmfdp,   zdpmel,   plu,  &
     &     plude,    ilab)

!----------------------------------
!*    3.0   cloud base calculations
!----------------------------------
!*         (a) determine cloud base values in 'cutypen',
!              and the cumulus type 1 or 2
!          -------------------------------------------
       call cutypen &
     &     ( klon,     klev,     klevp1,   klevm1,     pqen,&
     &      ztenh,    zqenh,     zqsenh,    zgeoh,     paph,&
     &      phhfl,    pqhfl,       pgeo,    pqsen,      pap,&
     &       pten,     lndj,        ptu,      pqu,     ilab,&
     &      ldcum,    kcbot,     ictop0,    ktype,    wbase,    plu,   kdpl)

!*         (b) assign the first guess mass flux at cloud base
!              ------------------------------------------
       do jl=1,klon
         zdhpbl(jl)=0.0
         upbl(jl) = 0.0
         idtop(jl)=0
       end do

       do jk=2,klev
       do jl=1,klon
         if(jk.ge.kcbot(jl) .and. ldcum(jl)) then
            zdhpbl(jl)=zdhpbl(jl)+(alv*pqte(jl,jk)+cpd*ptte(jl,jk))&
     &                 *(paph(jl,jk+1)-paph(jl,jk))
            if(lndj(jl) .eq. 0) then
              wspeed = sqrt(puen(jl,jk)**2 + pven(jl,jk)**2) 
              upbl(jl) = upbl(jl) + wspeed*(paph(jl,jk+1)-paph(jl,jk))
            end if
         end if
       end do
       end do

      do jl=1,klon
        if(ldcum(jl)) then
           ikb=kcbot(jl)
           zmfmax = (paph(jl,ikb)-paph(jl,ikb-1))*zcons2
           if(ktype(jl) == 1) then
             zmfub(jl)= 0.1*zmfmax
           else if ( ktype(jl) == 2 ) then
             zqumqe = pqu(jl,ikb) + plu(jl,ikb) - zqenh(jl,ikb)
             zdqmin = max(0.01*zqenh(jl,ikb),1.e-10)
             zdh = cpd*(ptu(jl,ikb)-ztenh(jl,ikb)) + alv*zqumqe
             zdh = g*max(zdh,1.e5*zdqmin)
             if ( zdhpbl(jl) > 0. ) then
               zmfub(jl) = zdhpbl(jl)/zdh
               zmfub(jl) = min(zmfub(jl),zmfmax)
             else
               zmfub(jl) = 0.1*zmfmax
               ldcum(jl) = .false.
             end if
            end if  
        else
           zmfub(jl) = 0.
        end if
      end do
!------------------------------------------------------
!*    4.0   determine cloud ascent for entraining plume
!------------------------------------------------------
!*    (a) do ascent in 'cuasc'in absence of downdrafts
!----------------------------------------------------------
      call cuascn &
     &    (klon,     klev,     klevp1,   klevm1,   ztenh,&
     &     zqenh,    puen,     pven,     pten,     pqen,&
     &     pqsen,    pgeo,     zgeoh,    pap,      paph,&
     &     pqte,     pverv,    ilwmin,   ldcum,    zhcbase,&
     &     ktype,    ilab,     ptu,      pqu,      plu,&
     &     zuu,      zvu,      pmfu,     zmfub,&
     &     zmfus,    zmfuq,    zmful,    plude,    zdmfup,&
     &     kcbot,    kctop,    ictop0,   icum,     ztmst,&
     &     zqsenh,   zlglac,   lndj,     wup,      wbase,   kdpl,  pmfude_rate )

!*     (b) check cloud depth and change entrainment rate accordingly
!          calculate precipitation rate (for downdraft calculation)
!------------------------------------------------------------------
      do jl=1,klon
        if ( ldcum(jl) ) then
          ikb = kcbot(jl)
          itopm2 = kctop(jl)
          zpbmpt = paph(jl,ikb) - paph(jl,itopm2)
          if ( ktype(jl) == 1 .and. zpbmpt <  zdnoprc ) ktype(jl) = 2
          if ( ktype(jl) == 2 .and. zpbmpt >= zdnoprc ) ktype(jl) = 1
          ictop0(jl) = kctop(jl)
        end if
        zrfl(jl)=zdmfup(jl,1)
      end do

      do jk=2,klev
        do jl=1,klon
          zrfl(jl)=zrfl(jl)+zdmfup(jl,jk)
        end do
      end do

      do jk = 1,klev
      do jl = 1,klon
        pmfd(jl,jk) = 0.
        zmfds(jl,jk) = 0.
        zmfdq(jl,jk) = 0.
        zdmfdp(jl,jk) = 0.
        zdpmel(jl,jk) = 0.
      end do
      end do

!-----------------------------------------
!*    6.0   cumulus downdraft calculations
!-----------------------------------------
      if(lmfdd) then
!*      (a) determine lfs in 'cudlfsn'
!--------------------------------------
        call cudlfsn    &
     &    (klon,     klev,&
     &     kcbot,    kctop,    lndj,   ldcum,   &
     &     ztenh,    zqenh,    puen,     pven,   &
     &     pten,     pqsen,    pgeo,              &
     &     zgeoh,    paph,     ptu,      pqu,      plu,   &
     &     zuu,      zvu,      zmfub,    zrfl,             &
     &     ztd,      zqd,      zud,      zvd,               &
     &     pmfd,     zmfds,    zmfdq,    zdmfdp,             &
     &     idtop,    loddraf)
!*     (b)  determine downdraft t,q and fluxes in 'cuddrafn'
!------------------------------------------------------------
        call cuddrafn                                         &
     &   ( klon,     klev,     loddraf,                      &
     &     ztenh,    zqenh,    puen,     pven,                 &
     &     pgeo,     zgeoh,    paph,     zrfl,                  &
     &     ztd,      zqd,      zud,      zvd,      pmfu,         &
     &     pmfd,     zmfds,    zmfdq,    zdmfdp,   pmfdde_rate      )
!-----------------------------------------------------------
      end if
!
!-----------------------------------------------------------------------
!* 6.0          closure and clean work
! ------
!-- 6.1 recalculate cloud base massflux from a cape closure
!       for deep convection (ktype=1) 
!
      do jl=1,klon
      if(ldcum(jl) .and. ktype(jl) .eq. 1) then
        ikb = kcbot(jl)
        ikt = kctop(jl)
        zheat(jl)=0.0
        zcape(jl)=0.0
        zcape1(jl)=0.0
        zcape2(jl)=0.0
        zmfub1(jl)=zmfub(jl)
    
        ztauc(jl)  = (zgeoh(jl,ikt)-zgeoh(jl,ikb)) / &
                   ((2.+ min(15.0,wup(jl)))*g)
        if(lndj(jl) .eq. 0) then 
          upbl(jl) = 2.+ upbl(jl)/(paph(jl,klev+1)-paph(jl,ikb))
          ztaubl(jl) = (zgeoh(jl,ikb)-zgeoh(jl,klev+1))/(g*upbl(jl))
          ztaubl(jl) = min(300., ztaubl(jl))
        else
          ztaubl(jl) = ztauc(jl)
        end if
      end if    
      end do
!
      do jk = 1 , klev
      do jl = 1 , klon
        llo1 = ldcum(jl) .and. ktype(jl) .eq. 1
        if ( llo1 .and. jk <= kcbot(jl) .and. jk > kctop(jl) ) then
          ikb = kcbot(jl)
          zdz = pgeo(jl,jk-1)-pgeo(jl,jk)
          zdp = pap(jl,jk)-pap(jl,jk-1)
          zheat(jl) = zheat(jl) + ((pten(jl,jk-1)-pten(jl,jk)+zdz*rcpd) / &
                      ztenh(jl,jk)+vtmpc1*(pqen(jl,jk-1)-pqen(jl,jk))) * &
                      (g*(pmfu(jl,jk)+pmfd(jl,jk)))
          zcape1(jl) = zcape1(jl) + ((ptu(jl,jk)-ztenh(jl,jk))/ztenh(jl,jk) + &
                      vtmpc1*(pqu(jl,jk)-zqenh(jl,jk))-plu(jl,jk))*zdp
        end if

        if ( llo1 .and. jk >= kcbot(jl) ) then
        if((paph(jl,klev+1)-paph(jl,kdpl(jl)))<50.e2) then
          zdp = paph(jl,jk+1)-paph(jl,jk)
          zcape2(jl) = zcape2(jl) + ztaubl(jl)* &
                     ((1.+vtmpc1*pqen(jl,jk))*ptte(jl,jk)+vtmpc1*pten(jl,jk)*pqte(jl,jk))*zdp 
        end if
        end if
      end do
      end do

      do jl=1,klon
       if(ldcum(jl).and.ktype(jl).eq.1) then
           ikb = kcbot(jl)
           ikt = kctop(jl)
           ztau = ztauc(jl) * (1.+1.33e-5*dx(jl))
           ztau = max(ztmst,ztau)
           ztau = max(720.,ztau)
           ztau = min(10800.,ztau)
           if(isequil) then
             zcape2(jl)= max(0.,zcape2(jl))
             zcape(jl) = max(0.,min(zcape1(jl)-zcape2(jl),5000.))
           else
             zcape(jl) = max(0.,min(zcape1(jl),5000.))
           end if
           zheat(jl) = max(1.e-4,zheat(jl))
           zmfub1(jl) = (zcape(jl)*zmfub(jl))/(zheat(jl)*ztau)
           zmfub1(jl) = max(zmfub1(jl),0.001)
           zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
           zmfub1(jl)=min(zmfub1(jl),zmfmax)
       end if
      end do
!
!*  6.2   recalculate convective fluxes due to effect of
!         downdrafts on boundary layer moist static energy budget (ktype=2)
!--------------------------------------------------------
       do jl=1,klon
         if(ldcum(jl) .and. ktype(jl) .eq. 2) then
           ikb=kcbot(jl)
           if(pmfd(jl,ikb).lt.0.0 .and. loddraf(jl)) then
              zeps=-pmfd(jl,ikb)/max(zmfub(jl),cmfcmin)
           else
              zeps=0.
           endif
           zqumqe=pqu(jl,ikb)+plu(jl,ikb)-  &
     &            zeps*zqd(jl,ikb)-(1.-zeps)*zqenh(jl,ikb)
           zdqmin=max(0.01*zqenh(jl,ikb),cmfcmin)
           zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
!  using moist static engergy closure instead of moisture closure
           zdh=cpd*(ptu(jl,ikb)-zeps*ztd(jl,ikb)- &
     &       (1.-zeps)*ztenh(jl,ikb))+alv*zqumqe
           zdh=g*max(zdh,1.e5*zdqmin)
           if(zdhpbl(jl).gt.0.)then
             zmfub1(jl)=zdhpbl(jl)/zdh
           else
             zmfub1(jl) = zmfub(jl)
           end if
           zmfub1(jl) = min(zmfub1(jl),zmfmax)
         end if

!*  6.3   mid-level convection - nothing special
!---------------------------------------------------------
         if(ldcum(jl) .and. ktype(jl) .eq. 3 ) then
            zmfub1(jl) = zmfub(jl)
         end if

       end do

!*  6.4   scaling the downdraft mass flux
!---------------------------------------------------------
       do jk=1,klev
       do jl=1,klon
        if( ldcum(jl) ) then
           zfac=zmfub1(jl)/max(zmfub(jl),cmfcmin)
           pmfd(jl,jk)=pmfd(jl,jk)*zfac
           zmfds(jl,jk)=zmfds(jl,jk)*zfac
           zmfdq(jl,jk)=zmfdq(jl,jk)*zfac
           zdmfdp(jl,jk)=zdmfdp(jl,jk)*zfac
           pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk)*zfac
        end if
       end do
       end do

!*  6.5   scaling the updraft mass flux
! --------------------------------------------------------
    do jl = 1,klon
      if ( ldcum(jl) ) zmfs(jl) = zmfub1(jl)/max(cmfcmin,zmfub(jl))
    end do
    do jk = 2 , klev
      do jl = 1,klon
        if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
          ikb = kcbot(jl)
          if ( jk>ikb ) then
            zdz = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ikb)))
            pmfu(jl,jk) = pmfu(jl,ikb)*zdz
          end if
          zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
          if ( pmfu(jl,jk)*zmfs(jl) > zmfmax ) then
            zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk))
          end if
        end if
      end do
    end do
    do jk = 2 , klev
      do jl = 1,klon
        if ( ldcum(jl) .and. jk <= kcbot(jl) .and. jk >= kctop(jl)-1 ) then
          pmfu(jl,jk) = pmfu(jl,jk)*zmfs(jl)
          zmfus(jl,jk) = zmfus(jl,jk)*zmfs(jl)
          zmfuq(jl,jk) = zmfuq(jl,jk)*zmfs(jl)
          zmful(jl,jk) = zmful(jl,jk)*zmfs(jl)
          zdmfup(jl,jk) = zdmfup(jl,jk)*zmfs(jl)
          plude(jl,jk) = plude(jl,jk)*zmfs(jl)
          pmfude_rate(jl,jk) = pmfude_rate(jl,jk)*zmfs(jl)
        end if
      end do
    end do

!*    6.6  if ktype = 2, kcbot=kctop is not allowed
! ---------------------------------------------------
    do jl = 1,klon
      if ( ktype(jl) == 2 .and. &
           kcbot(jl) == kctop(jl) .and. kcbot(jl) >= klev-1 ) then
        ldcum(jl) = .false.
        ktype(jl) = 0
      end if
    end do

    if ( .not. lmfscv .or. .not. lmfpen ) then
      do jl = 1,klon
        llo2(jl) = .false.
        if ( (.not. lmfscv .and. ktype(jl) == 2) .or. &
             (.not. lmfpen .and. ktype(jl) == 1) ) then
          llo2(jl) = .true.
          ldcum(jl) = .false.
        end if
      end do
    end if

!*   6.7  set downdraft mass fluxes to zero above cloud top
!----------------------------------------------------
    do jl = 1,klon
      if ( loddraf(jl) .and. idtop(jl) <= kctop(jl) ) then
        idtop(jl) = kctop(jl) + 1
      end if
    end do
    do jk = 2 , klev
      do jl = 1,klon
        if ( loddraf(jl) ) then
          if ( jk < idtop(jl) ) then
            pmfd(jl,jk) = 0.
            zmfds(jl,jk) = 0.
            zmfdq(jl,jk) = 0.
            pmfdde_rate(jl,jk) = 0.
            zdmfdp(jl,jk) = 0.
          else if ( jk == idtop(jl) ) then
            pmfdde_rate(jl,jk) = 0.
          end if
        end if
      end do
    end do

    itopm2 = 2
!----------------------------------------------------------
!*    7.0      determine final convective fluxes in 'cuflx'
!----------------------------------------------------------
       call cuflxn                                        &                  
     &  (  klon,     klev,     ztmst                       &                                       
     &  ,  pten,     pqen,     pqsen,    ztenh,   zqenh     &                  
     &  ,  paph,     pap,      zgeoh,    lndj,    ldcum      &                 
     &  ,  kcbot,    kctop,    idtop,    itopm2               &                 
     &  ,  ktype,    loddraf                                   &                 
     &  ,  pmfu,     pmfd,     zmfus,    zmfds                  &               
     &  ,  zmfuq,    zmfdq,    zmful,    plude                   &              
     &  ,  zdmfup,   zdmfdp,   zdpmel,   zlglac                   &             
     &  ,  prain,    pmfdde_rate, pmflxr, pmflxs )     

! some adjustments needed
    do jl=1,klon
      zmfs(jl) = 1.
      zmfuub(jl)=0.
    end do
    do jk = 2 , klev
      do jl = 1,klon
        if ( loddraf(jl) .and. jk >= idtop(jl)-1 ) then
          zmfmax = pmfu(jl,jk)*0.98
          if ( pmfd(jl,jk)+zmfmax+1.e-15 < 0. ) then
            zmfs(jl) = min(zmfs(jl),-zmfmax/pmfd(jl,jk))
          end if
        end if
      end do
    end do

    do jk = 2 , klev
      do jl = 1 , klon
        if ( zmfs(jl) < 1. .and. jk >= idtop(jl)-1 ) then
          pmfd(jl,jk) = pmfd(jl,jk)*zmfs(jl)
          zmfds(jl,jk) = zmfds(jl,jk)*zmfs(jl)
          zmfdq(jl,jk) = zmfdq(jl,jk)*zmfs(jl)
          pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk)*zmfs(jl)
          zmfuub(jl) = zmfuub(jl) - (1.-zmfs(jl))*zdmfdp(jl,jk)
          pmflxr(jl,jk+1) = pmflxr(jl,jk+1) + zmfuub(jl)
          zdmfdp(jl,jk) = zdmfdp(jl,jk)*zmfs(jl)
        end if
      end do
    end do

    do jk = 2 , klev - 1
      do jl = 1, klon
        if ( loddraf(jl) .and. jk >= idtop(jl)-1 ) then
          zerate = -pmfd(jl,jk) + pmfd(jl,jk-1) + pmfdde_rate(jl,jk)
          if ( zerate < 0. ) then
            pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk) - zerate
          end if
        end if
        if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
          zerate = pmfu(jl,jk) - pmfu(jl,jk+1) + pmfude_rate(jl,jk)
          if ( zerate < 0. ) then
            pmfude_rate(jl,jk) = pmfude_rate(jl,jk) - zerate
          end if
          zdmfup(jl,jk) = pmflxr(jl,jk+1) + pmflxs(jl,jk+1) - &
                          pmflxr(jl,jk) - pmflxs(jl,jk)
          zdmfdp(jl,jk) = 0.
        end if
      end do
    end do

! avoid negative humidities at ddraught top
    do jl = 1,klon
      if ( loddraf(jl) ) then
        jk = idtop(jl)
        ik = min(jk+1,klev)
        if ( zmfdq(jl,jk) < 0.3*zmfdq(jl,ik) ) then
            zmfdq(jl,jk) = 0.3*zmfdq(jl,ik)
        end if
      end if
    end do

! avoid negative humidities near cloud top because gradient of precip flux
! and detrainment / liquid water flux are too large
    do jk = 2 , klev
      do jl = 1, klon
        if ( ldcum(jl) .and. jk >= kctop(jl)-1 .and. jk < kcbot(jl) ) then
          zdz = ztmst*g/(paph(jl,jk+1)-paph(jl,jk))
          zmfa = zmfuq(jl,jk+1) + zmfdq(jl,jk+1) - &
                 zmfuq(jl,jk) - zmfdq(jl,jk) + &
                 zmful(jl,jk+1) - zmful(jl,jk) + zdmfup(jl,jk)
          zmfa = (zmfa-plude(jl,jk))*zdz
          if ( pqen(jl,jk)+zmfa < 0. ) then
            plude(jl,jk) = plude(jl,jk) + 2.*(pqen(jl,jk)+zmfa)/zdz
          end if
          if ( plude(jl,jk) < 0. ) plude(jl,jk) = 0.
        end if
        if ( .not. ldcum(jl) ) pmfude_rate(jl,jk) = 0.
        if ( abs(pmfd(jl,jk-1)) < 1.0e-20 ) pmfdde_rate(jl,jk) = 0.
      end do
    end do

    do jl=1,klon
      prsfc(jl) = pmflxr(jl,klev+1)
      pssfc(jl) = pmflxs(jl,klev+1)
    end do

!----------------------------------------------------------------
!*    8.0      update tendencies for t and q in subroutine cudtdq
!----------------------------------------------------------------
      call cudtdqn(klon,klev,itopm2,kctop,idtop,ldcum,loddraf, &
                 ztmst,paph,zgeoh,pgeo,pten,ztenh,pqen,zqenh,pqsen,     &
                 zlglac,plude,pmfu,pmfd,zmfus,zmfds,zmfuq,zmfdq,zmful,   &
                 zdmfup,zdmfdp,zdpmel,ptte,pqte,pcte)
!----------------------------------------------------------------
!*    9.0      update tendencies for u and u in subroutine cududv
!----------------------------------------------------------------
      if(lmfdudv) then
      do jk = klev-1 , 2 , -1
        ik = jk + 1
        do jl = 1,klon
          if ( ldcum(jl) ) then
            if ( jk == kcbot(jl) .and. ktype(jl) < 3 ) then
              ikb = kdpl(jl)
              zuu(jl,jk) = puen(jl,ikb-1)
              zvu(jl,jk) = pven(jl,ikb-1)
            else if ( jk == kcbot(jl) .and. ktype(jl) == 3 ) then
              zuu(jl,jk) = puen(jl,jk-1)
              zvu(jl,jk) = pven(jl,jk-1)
            end if
            if ( jk < kcbot(jl) .and. jk >= kctop(jl) ) then
            if(momtrans .eq. 1)then
              zfac = 0.
              if ( ktype(jl) == 1 .or. ktype(jl) == 3 ) zfac = 2.
              if ( ktype(jl) == 1 .and. jk <= kctop(jl)+2 ) zfac = 3.
              zerate = pmfu(jl,jk) - pmfu(jl,ik) + &
                (1.+zfac)*pmfude_rate(jl,jk)
              zderate = (1.+zfac)*pmfude_rate(jl,jk)
              zmfa = 1./max(cmfcmin,pmfu(jl,jk))
              zuu(jl,jk) = (zuu(jl,ik)*pmfu(jl,ik) + &
                zerate*puen(jl,jk)-zderate*zuu(jl,ik))*zmfa
              zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + &
                zerate*pven(jl,jk)-zderate*zvu(jl,ik))*zmfa
            else
              if(ktype(jl) == 1 .or. ktype(jl) == 3) then
                pgf_u   =  -0.7*0.5*(pmfu(jl,ik)*(puen(jl,ik)-puen(jl,jk))+&
                                   pmfu(jl,jk)*(puen(jl,jk)-puen(jl,jk-1)))
                pgf_v   =  -0.7*0.5*(pmfu(jl,ik)*(pven(jl,ik)-pven(jl,jk))+&
                                   pmfu(jl,jk)*(pven(jl,jk)-pven(jl,jk-1)))
              else
                pgf_u   = 0.
                pgf_v   = 0.
              end if
              zerate = pmfu(jl,jk) - pmfu(jl,ik) + pmfude_rate(jl,jk)
              zderate = pmfude_rate(jl,jk)
              zmfa = 1./max(cmfcmin,pmfu(jl,jk))
              zuu(jl,jk) = (zuu(jl,ik)*pmfu(jl,ik) + &
                zerate*puen(jl,jk)-zderate*zuu(jl,ik)+pgf_u)*zmfa
              zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + &
                zerate*pven(jl,jk)-zderate*zvu(jl,ik)+pgf_v)*zmfa
            end if
            end if
          end if
        end do
      end do

      if(lmfdd) then
      do jk = 3 , klev
        ik = jk - 1
        do jl = 1,klon
          if ( ldcum(jl) ) then
            if ( jk == idtop(jl) ) then
              zud(jl,jk) = 0.5*(zuu(jl,jk)+puen(jl,ik))
              zvd(jl,jk) = 0.5*(zvu(jl,jk)+pven(jl,ik))
            else if ( jk > idtop(jl) ) then
              zerate = -pmfd(jl,jk) + pmfd(jl,ik) + pmfdde_rate(jl,jk)
              zmfa = 1./min(-cmfcmin,pmfd(jl,jk))
              zud(jl,jk) = (zud(jl,ik)*pmfd(jl,ik) - &
                zerate*puen(jl,ik)+pmfdde_rate(jl,jk)*zud(jl,ik))*zmfa
              zvd(jl,jk) = (zvd(jl,ik)*pmfd(jl,ik) - &
                zerate*pven(jl,ik)+pmfdde_rate(jl,jk)*zvd(jl,ik))*zmfa
            end if
          end if
        end do
      end do
      end if
!   --------------------------------------------------
!   rescale massfluxes for stability in Momentum
!------------------------------------------------------------------------
      zmfs(:) = 1.
      do jk = 2 , klev
        do jl = 1, klon
          if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
            zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons
            if ( pmfu(jl,jk) > zmfmax .and. jk >= kctop(jl) ) then
              zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk))
            end if
          end if
        end do
      end do
      do jk = 1 , klev
        do jl = 1, klon
          zmfuus(jl,jk) = pmfu(jl,jk)
          zmfdus(jl,jk) = pmfd(jl,jk)
          if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
            zmfuus(jl,jk) = pmfu(jl,jk)*zmfs(jl)
            zmfdus(jl,jk) = pmfd(jl,jk)*zmfs(jl)
          end if
        end do
      end do
!*  9.1          update u and v in subroutine cududvn
!-------------------------------------------------------------------
     do jk = 1 , klev
        do jl = 1, klon
          ztenu(jl,jk) = pvom(jl,jk)
          ztenv(jl,jk) = pvol(jl,jk)
        end do
      end do

      call cududvn(klon,klev,itopm2,ktype,kcbot,kctop, &
                  ldcum,ztmst,paph,puen,pven,zmfuus,zmfdus,zuu,  &
                  zud,zvu,zvd,pvom,pvol)

!  calculate KE dissipation
      do jl = 1, klon
        zsum12(jl) = 0.
        zsum22(jl) = 0.
      end do
        do jk = 1 , klev
          do jl = 1, klon
            zuv2(jl,jk) = 0.
            if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
              zdz = (paph(jl,jk+1)-paph(jl,jk))
              zduten = pvom(jl,jk) - ztenu(jl,jk)
              zdvten = pvol(jl,jk) - ztenv(jl,jk)
              zuv2(jl,jk) = sqrt(zduten**2+zdvten**2)
              zsum22(jl) = zsum22(jl) + zuv2(jl,jk)*zdz
              zsum12(jl) = zsum12(jl) - &
                (puen(jl,jk)*zduten+pven(jl,jk)*zdvten)*zdz
            end if
          end do
        end do
        do jk = 1 , klev
          do jl = 1, klon
            if ( ldcum(jl) .and. jk>=kctop(jl)-1 ) then
              ztdis = rcpd*zsum12(jl)*zuv2(jl,jk)/max(1.e-15,zsum22(jl))
              ptte(jl,jk) = ptte(jl,jk) + ztdis
            end if
          end do
        end do

      end if

!----------------------------------------------------------------------
!*   10.           IN CASE THAT EITHER DEEP OR SHALLOW IS SWITCHED OFF
! NEED TO SET SOME VARIABLES A POSTERIORI TO ZERO
! ---------------------------------------------------
    if ( .not. lmfscv .or. .not. lmfpen ) then
      do jk = 2 , klev
        do jl = 1, klon
          if ( llo2(jl) .and. jk >= kctop(jl)-1 ) then
            ptu(jl,jk) = pten(jl,jk)
            pqu(jl,jk) = pqen(jl,jk)
            plu(jl,jk) = 0.
            pmfude_rate(jl,jk) = 0.
            pmfdde_rate(jl,jk) = 0.
          end if
        end do
      end do
      do jl = 1, klon
        if ( llo2(jl) ) then
          kctop(jl) = klev - 1
          kcbot(jl) = klev - 1
        end if
      end do
    end if

 !----------------------------------------------------------------------
 !*   11.0          CHEMICAL TRACER TRANSPORT
 ! -------------------------

    if (  ktrac > 0 ) then
      !   transport switched off for mid-level convection
      do jl = 1, klon
        if ( ldcum(jl) .and. ktype(jl) /= 3 .and. &
             kcbot(jl)-kctop(jl) >= 1 ) then
          lldcum(jl) = .true.
          llddraf3(jl) = loddraf(jl)
        else
          lldcum(jl) = .false.
          llddraf3(jl) = .false.
        end if
      end do
      !   check and correct mass fluxes for CFL criterium
      zmfs(:) = 1.
      do jk = 2 , klev
        do jl = 1, klon
          if ( lldcum(jl) .and. jk >= kctop(jl) ) then
            zmfmax = (paph(jl,jk)-paph(jl,jk-1))*0.8*zcons
            if ( pmfu(jl,jk) > zmfmax ) then
              zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk))
            end if
          end if
        end do
      end do

      do jk = 1, klev
        do jl = 1, klon
          if ( lldcum(jl) .and. jk >= kctop(jl)-1 ) then
            zmfuus(jl,jk) = pmfu(jl,jk)*zmfs(jl)
            zmfudr(jl,jk) = pmfude_rate(jl,jk)*zmfs(jl)
          else
            zmfuus(jl,jk) = 0.
            zmfudr(jl,jk) = 0.
          end if
          if ( llddraf3(jl) .and. jk >= idtop(jl)-1 ) then
            zmfdus(jl,jk) = pmfd(jl,jk)*zmfs(jl)
            zmfddr(jl,jk) = pmfdde_rate(jl,jk)*zmfs(jl)
          else
            zmfdus(jl,jk) = 0.
            zmfddr(jl,jk) = 0.
          end if
        end do
      end do

      call cuctracer(klon,klev,ktrac,kctop,idtop,   &
                     lldcum,llddraf3,ztmst,paph,zmfuus,zmfdus, &
                     zmfudr,zmfddr,pcen,ptenc)
    end if

      return
      end subroutine cumastrn

!**********************************************
!    level 3  subroutine cuinin
!**********************************************
!
      subroutine cuinin &
     &    (klon,     klev,     klevp1,   klevm1,   pten,&
     &     pqen,     pqsen,    puen,     pven,     pverv,&
     &     pgeo,     paph,     pgeoh,    ptenh,    pqenh,&
     &     pqsenh,   klwmin,   ptu,      pqu,      ptd,&
     &     pqd,      puu,      pvu,      pud,      pvd,&
     &     pmfu,     pmfd,     pmfus,    pmfds,    pmfuq,&
     &     pmfdq,    pdmfup,   pdmfdp,   pdpmel,   plu,&
     &     plude,    klab)
      implicit none
!      m.tiedtke         e.c.m.w.f.     12/89
!***purpose
!   -------
!          this routine interpolates large-scale fields of t,q etc.
!          to half levels (i.e. grid for massflux scheme),
!          and initializes values for updrafts and downdrafts
!***interface
!   ---------
!          this routine is called from *cumastr*.
!***method.
!  --------
!          for extrapolation to half levels see tiedtke(1989)
!***externals
!   ---------
!          *cuadjtq* to specify qs at half levels
! ----------------------------------------------------------------
      integer  klon,klev,klevp1,klevm1
      real(kind=kind_phys)     pten(klon,klev),        pqen(klon,klev),&
     &         puen(klon,klev),        pven(klon,klev),&
     &         pqsen(klon,klev),       pverv(klon,klev),&
     &         pgeo(klon,klev),        pgeoh(klon,klevp1),&
     &         paph(klon,klevp1),      ptenh(klon,klev),&
     &         pqenh(klon,klev),       pqsenh(klon,klev)
      real(kind=kind_phys)     ptu(klon,klev),         pqu(klon,klev),&
     &         ptd(klon,klev),         pqd(klon,klev),&
     &         puu(klon,klev),         pud(klon,klev),&
     &         pvu(klon,klev),         pvd(klon,klev),&
     &         pmfu(klon,klev),        pmfd(klon,klev),&
     &         pmfus(klon,klev),       pmfds(klon,klev),&
     &         pmfuq(klon,klev),       pmfdq(klon,klev),&
     &         pdmfup(klon,klev),      pdmfdp(klon,klev),&
     &         plu(klon,klev),         plude(klon,klev)
      real(kind=kind_phys)     zwmax(klon),            zph(klon), &
     &         pdpmel(klon,klev)
      integer  klab(klon,klev),        klwmin(klon)
      logical  loflag(klon)
!  local variables
      integer  jl,jk
      integer  icall,ik
      real(kind=kind_phys)     zzs
!------------------------------------------------------------
!*    1.       specify large scale parameters at half levels
!*             adjust temperature fields if staticly unstable
!*             find level of maximum vertical velocity
! -----------------------------------------------------------
      do jk=2,klev
      do jl=1,klon
        ptenh(jl,jk)=(max(cpd*pten(jl,jk-1)+pgeo(jl,jk-1), &
     &             cpd*pten(jl,jk)+pgeo(jl,jk))-pgeoh(jl,jk))*rcpd
        pqenh(jl,jk) = pqen(jl,jk-1)
        pqsenh(jl,jk)= pqsen(jl,jk-1)
        zph(jl)=paph(jl,jk)
        loflag(jl)=.true.
      end do

      if ( jk >= klev-1 .or. jk < 2 ) cycle
      ik=jk
      icall=0
      call cuadjtqn(klon,klev,ik,zph,ptenh,pqsenh,loflag,icall)
      do jl=1,klon
        pqenh(jl,jk)=min(pqen(jl,jk-1),pqsen(jl,jk-1)) &
     &            +(pqsenh(jl,jk)-pqsen(jl,jk-1))
        pqenh(jl,jk)=max(pqenh(jl,jk),0.)
      end do
      end do

      do jl=1,klon
        ptenh(jl,klev)=(cpd*pten(jl,klev)+pgeo(jl,klev)- &
     &                pgeoh(jl,klev))*rcpd
        pqenh(jl,klev)=pqen(jl,klev)
        ptenh(jl,1)=pten(jl,1)
        pqenh(jl,1)=pqen(jl,1)
        klwmin(jl)=klev
        zwmax(jl)=0.
      end do

      do jk=klevm1,2,-1
      do jl=1,klon
        zzs=max(cpd*ptenh(jl,jk)+pgeoh(jl,jk), &
     &        cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))
        ptenh(jl,jk)=(zzs-pgeoh(jl,jk))*rcpd
      end do
      end do

      do jk=klev,3,-1
      do jl=1,klon
        if(pverv(jl,jk).lt.zwmax(jl)) then
           zwmax(jl)=pverv(jl,jk)
           klwmin(jl)=jk
        end if
      end do
      end do
!-----------------------------------------------------------
!*    2.0      initialize values for updrafts and downdrafts
!-----------------------------------------------------------
      do jk=1,klev
      ik=jk-1
      if(jk.eq.1) ik=1
      do jl=1,klon
      ptu(jl,jk)=ptenh(jl,jk)
      ptd(jl,jk)=ptenh(jl,jk)
      pqu(jl,jk)=pqenh(jl,jk)
      pqd(jl,jk)=pqenh(jl,jk)
      plu(jl,jk)=0.
      puu(jl,jk)=puen(jl,ik)
      pud(jl,jk)=puen(jl,ik)
      pvu(jl,jk)=pven(jl,ik)
      pvd(jl,jk)=pven(jl,ik)
      klab(jl,jk)=0
      end do
      end do
      return
      end subroutine cuinin

!---------------------------------------------------------
!  level 3 souroutines
!--------------------------------------------------------
      subroutine cutypen & 
     &   (  klon,    klev,     klevp1,   klevm1,     pqen,&
     &     ptenh,   pqenh,     pqsenh,    pgeoh,     paph,& 
     &       hfx,     qfx,       pgeo,    pqsen,      pap,&
     &      pten,    lndj,       cutu,     cuqu,    culab,&
     &     ldcum,   cubot,      cutop,    ktype,    wbase,    culu,   kdpl  )
!      zhang & wang      iprc           2011-2013
!***purpose.
!   --------
!          to produce first guess updraught for cu-parameterizations
!          calculates condensation level, and sets updraught base variables and
!          first guess cloud type
!***interface
!   ---------
!          this routine is called from *cumastr*.
!          input are environm. values of t,q,p,phi at half levels.
!          it returns cloud types as follows;
!                 ktype=1 for deep cumulus
!                 ktype=2 for shallow cumulus
!***method.
!  --------
!          based on a simplified updraught equation
!            partial(hup)/partial(z)=eta(h - hup)
!            eta is the entrainment rate for test parcel
!            h stands for dry static energy or the total water specific humidity
!            references: christian jakob, 2003: a new subcloud model for
!            mass-flux convection schemes
!                        influence on triggering, updraft properties, and model
!                        climate, mon.wea.rev.
!                        131, 2765-2778
!            and
!                        ifs documentation - cy36r1,cy38r1 
!***input variables:
!       ptenh [ztenh] - environment temperature on half levels. (cuini)
!       pqenh [zqenh] - env. specific humidity on half levels. (cuini)
!       pgeoh [zgeoh] - geopotential on half levels, (mssflx)
!       paph - pressure of half levels. (mssflx)
!       rho  - density of the lowest model level
!       qfx  - net upward moisture flux at the surface (kg/m^2/s)
!       hfx  - net upward heat flux at the surface (w/m^2)
!***variables output by cutype:
!       ktype - convection type - 1: penetrative  (cumastr)
!                                 2: stratocumulus (cumastr)
!                                 3: mid-level  (cuasc)
!       information for updraft parcel (ptu,pqu,plu,kcbot,klab,kdpl...)
! ----------------------------------------------------------------
!-------------------------------------------------------------------
      implicit none
!-------------------------------------------------------------------
      integer  klon, klev, klevp1, klevm1
      real(kind=kind_phys)     ptenh(klon,klev),       pqenh(klon,klev),& 
     &         pqsen(klon,klev),       pqsenh(klon,klev),&
     &         pgeoh(klon,klevp1),     paph(klon,klevp1),&
     &         pap(klon,klev),         pqen(klon,klev)
      real(kind=kind_phys)     pten(klon,klev)
      real(kind=kind_phys)     ptu(klon,klev),pqu(klon,klev),plu(klon,klev)
      real(kind=kind_phys)     pgeo(klon,klev)
      integer  klab(klon,klev)
      integer  kctop(klon),kcbot(klon)

      real(kind=kind_phys)     qfx(klon),hfx(klon)
      real(kind=kind_phys)     zph(klon)
      integer  lndj(klon)
      logical  loflag(klon), deepflag(klon), resetflag(klon)

! output variables
      real(kind=kind_phys)     cutu(klon,klev), cuqu(klon,klev), culu(klon,klev)
      integer  culab(klon,klev)
      real(kind=kind_phys)     wbase(klon)
      integer  ktype(klon),cubot(klon),cutop(klon),kdpl(klon)
      logical  ldcum(klon)

! local variables
      real(kind=kind_phys)     zqold(klon)
      real(kind=kind_phys)     rho, part1, part2, root, conw, deltt, deltq
      real(kind=kind_phys)     eta(klon),dz(klon),coef(klon)
      real(kind=kind_phys)     dhen(klon,klev), dh(klon,klev)
      real(kind=kind_phys)     plude(klon,klev)
      real(kind=kind_phys)     kup(klon,klev)
      real(kind=kind_phys)     vptu(klon,klev),vten(klon,klev)
      real(kind=kind_phys)     zbuo(klon,klev),abuoy(klon,klev)
  
      real(kind=kind_phys)     zz,zdken,zdq
      real(kind=kind_phys)     fscale,crirh1,pp
      real(kind=kind_phys)     atop1,atop2,abot
      real(kind=kind_phys)     tmix,zmix,qmix,pmix
      real(kind=kind_phys)     zlglac,dp
      integer  nk,is,ikb,ikt

      real(kind=kind_phys)     zqsu,zcor,zdp,zesdp,zalfaw,zfacw,zfaci,zfac,zdsdp,zdqsdt,zdtdp
      real(kind=kind_phys)     zpdifftop, zpdiffbot
      integer  zcbase(klon), itoppacel(klon)
      integer  jl,jk,ik,icall,levels
      logical  needreset, lldcum(klon)
!--------------------------------------------------------------
      do jl=1,klon
        kcbot(jl)=klev
        kctop(jl)=klev
        kdpl(jl) =klev
        ktype(jl)=0
        wbase(jl)=0.
        ldcum(jl)=.false.
      end do

!-----------------------------------------------------------
! let's do test,and check the shallow convection first
! the first level is klev
! define deltat and deltaq
!-----------------------------------------------------------
      do jk=1,klev
      do jl=1,klon
          plu(jl,jk)=culu(jl,jk)  ! parcel liquid water
          ptu(jl,jk)=cutu(jl,jk)  ! parcel temperature
          pqu(jl,jk)=cuqu(jl,jk)  ! parcel specific humidity
          klab(jl,jk)=culab(jl,jk)
           dh(jl,jk)=0.0  ! parcel dry static energy
         dhen(jl,jk)=0.0  ! environment dry static energy
          kup(jl,jk)=0.0  ! updraught kinetic energy for parcel
         vptu(jl,jk)=0.0  ! parcel virtual temperature considering water-loading
         vten(jl,jk)=0.0  ! environment virtual temperature
         zbuo(jl,jk)=0.0  ! parcel buoyancy
         abuoy(jl,jk)=0.0
      end do
      end do

      do jl=1,klon
         zqold(jl) = 0.
         lldcum(jl) = .false.
         loflag(jl) = .true.
      end do

! check the levels from lowest level to second top level
      do jk=klevm1,2,-1

! define the variables at the first level      
      if(jk .eq. klevm1) then
      do jl=1,klon
        rho=pap(jl,klev)/ &
     &         (rd*(pten(jl,klev)*(1.+vtmpc1*pqen(jl,klev))))
        hfx(jl) = hfx(jl)*rho*cpd
        qfx(jl) = qfx(jl)*rho
        part1 = 1.5*0.4*pgeo(jl,klev)/ &
     &              (rho*pten(jl,klev))
        part2 = -hfx(jl)*rcpd-vtmpc1*pten(jl,klev)*qfx(jl)
        root  = 0.001-part1*part2
        if(part2 .lt. 0.) then
           conw  = 1.2*(root)**t13
           deltt = max(1.5*hfx(jl)/(rho*cpd*conw),0.)
           deltq = max(1.5*qfx(jl)/(rho*conw),0.)
           kup(jl,klev) = 0.5*(conw**2)
           pqu(jl,klev)= pqenh(jl,klev) + deltq
          dhen(jl,klev)= pgeoh(jl,klev) + ptenh(jl,klev)*cpd
           dh(jl,klev) = dhen(jl,klev)  + deltt*cpd
          ptu(jl,klev) = (dh(jl,klev)-pgeoh(jl,klev))*rcpd 
          vptu(jl,klev)=ptu(jl,klev)*(1.+vtmpc1*pqu(jl,klev))
          vten(jl,klev)=ptenh(jl,klev)*(1.+vtmpc1*pqenh(jl,klev))
          zbuo(jl,klev)=(vptu(jl,klev)-vten(jl,klev))/vten(jl,klev)
          klab(jl,klev) = 1
        else
          loflag(jl) = .false.
        end if
      end do
      end if
 
      is=0
      do jl=1,klon
         if(loflag(jl))then
            is=is+1
         endif
      enddo
      if(is.eq.0) exit

! the next levels, we use the variables at the first level as initial values
      do jl=1,klon
      if(loflag(jl)) then
        eta(jl) = 0.55/(pgeo(jl,jk)*zrg)+1.e-4
        dz(jl)  = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
        coef(jl)= 0.5*eta(jl)*dz(jl)
        dhen(jl,jk) = pgeoh(jl,jk) + cpd*ptenh(jl,jk)
        dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk))&
     &              +(1.-coef(jl))*dh(jl,jk+1))/(1.+coef(jl))
        pqu(jl,jk) =(coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk))&
     &              +(1.-coef(jl))*pqu(jl,jk+1))/(1.+coef(jl))
        ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd
        zqold(jl) = pqu(jl,jk)
        zph(jl)=paph(jl,jk)
      end if
      end do
! check if the parcel is saturated
      ik=jk
      icall=1
      call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
      do jl=1,klon
        if( loflag(jl) ) then
          zdq = max((zqold(jl) - pqu(jl,jk)),0.)
          plu(jl,jk) = plu(jl,jk+1) + zdq
          zlglac=zdq*((1.-foealfa(ptu(jl,jk))) - &
                      (1.-foealfa(ptu(jl,jk+1))))
          plu(jl,jk) = min(plu(jl,jk),5.e-3)
          dh(jl,jk) =  pgeoh(jl,jk) + cpd*(ptu(jl,jk)+ralfdcp*zlglac)
! compute the updraft speed
          vptu(jl,jk) = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))+&
                        ralfdcp*zlglac
          vten(jl,jk) = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
          zbuo(jl,jk) = (vptu(jl,jk) - vten(jl,jk))/vten(jl,jk)
          abuoy(jl,jk)=(zbuo(jl,jk)+zbuo(jl,jk+1))*0.5*g
          atop1 = 1.0 - 2.*coef(jl)
          atop2 = 2.0*dz(jl)*abuoy(jl,jk)
          abot =  1.0 + 2.*coef(jl)
          kup(jl,jk)  = (atop1*kup(jl,jk+1) + atop2) / abot

! let's find the exact cloud base
         if ( plu(jl,jk) > 0. .and. klab(jl,jk+1) == 1 ) then
              ik = jk + 1
              zqsu = foeewm(ptu(jl,ik))/paph(jl,ik)
              zqsu = min(0.5,zqsu)
              zcor = 1./(1.-vtmpc1*zqsu)
              zqsu = zqsu*zcor
              zdq = min(0.,pqu(jl,ik)-zqsu)
              zalfaw = foealfa(ptu(jl,ik))
              zfacw = c5les/((ptu(jl,ik)-c4les)**2)
              zfaci = c5ies/((ptu(jl,ik)-c4ies)**2)
              zfac = zalfaw*zfacw + (1.-zalfaw)*zfaci
              zesdp = foeewm(ptu(jl,ik))/paph(jl,ik)
              zcor = 1./(1.-vtmpc1*zesdp)
              zdqsdt = zfac*zcor*zqsu
              zdtdp = rd*ptu(jl,ik)/(cpd*paph(jl,ik))
              zdp = zdq/(zdqsdt*zdtdp)
              zcbase(jl) = paph(jl,ik) + zdp
! chose nearest half level as cloud base (jk or jk+1)
              zpdifftop = zcbase(jl) - paph(jl,jk)
              zpdiffbot = paph(jl,jk+1) - zcbase(jl)
              if ( zpdifftop > zpdiffbot .and. kup(jl,jk+1) > 0. ) then
                ikb = min(klev-1,jk+1)
                klab(jl,ikb) = 2
                klab(jl,jk) = 2
                kcbot(jl) = ikb
                plu(jl,jk+1) = 1.0e-8
              else if ( zpdifftop <= zpdiffbot .and.kup(jl,jk) > 0. ) then
                klab(jl,jk) = 2
                kcbot(jl) = jk
              end if
          end if

          if(kup(jl,jk) .lt. 0.)then
            loflag(jl) = .false.
            if(plu(jl,jk+1) .gt. 0.) then
              kctop(jl) = jk
              lldcum(jl) = .true.
            else
              lldcum(jl) = .false.
            end if
          else 
            if(plu(jl,jk) .gt. 0.)then
              klab(jl,jk)=2
            else
              klab(jl,jk)=1
            end if
          end if
        end if
      end do

      end do ! end all the levels

      do jl=1,klon
        ikb = kcbot(jl)
        ikt = kctop(jl)
        if(paph(jl,ikb) - paph(jl,ikt) > zdnoprc) lldcum(jl) = .false.
        if(lldcum(jl)) then
            ktype(jl) = 2
            ldcum(jl) = .true.
            wbase(jl) = sqrt(max(2.*kup(jl,ikb),0.))
            cubot(jl) = ikb
            cutop(jl) = ikt
            kdpl(jl)  = klev
        else
            cutop(jl) = -1
            cubot(jl) = -1
            kdpl(jl)  = klev - 1
            ldcum(jl) = .false.
            wbase(jl) = 0.
        end if
      end do

      do jk=klev,1,-1
       do jl=1,klon
           ikt = kctop(jl)
           if(jk .ge. ikt)then
             culab(jl,jk) = klab(jl,jk)
             cutu(jl,jk)  = ptu(jl,jk)
             cuqu(jl,jk)  = pqu(jl,jk)
             culu(jl,jk)  = plu(jl,jk)
           end if
       end do
      end do
      
!-----------------------------------------------------------
! next, let's check the deep convection
! the first level is klevm1-1
! define deltat and deltaq
!----------------------------------------------------------
! we check the parcel starting level by level
! assume the mix-layer is 60hPa
      deltt = 0.2
      deltq = 1.0e-4
      do jl=1,klon
        deepflag(jl) = .false.
      end do

      do jk=klev,1,-1
       do jl=1,klon
         if((paph(jl,klev+1)-paph(jl,jk)) .lt. 350.e2) itoppacel(jl) = jk
       end do
      end do

      do levels=klevm1-1,klevm1-20,-1 ! loop starts
        do jk=1,klev
          do jl=1,klon
             plu(jl,jk)=0.0  ! parcel liquid water
             ptu(jl,jk)=0.0  ! parcel temperature
             pqu(jl,jk)=0.0  ! parcel specific humidity
             dh(jl,jk)=0.0   ! parcel dry static energy
             dhen(jl,jk)=0.0  ! environment dry static energy
             kup(jl,jk)=0.0   ! updraught kinetic energy for parcel
             vptu(jl,jk)=0.0  ! parcel virtual temperature consideringwater-loading
             vten(jl,jk)=0.0  ! environment virtual temperature
             abuoy(jl,jk)=0.0
             zbuo(jl,jk)=0.0
             klab(jl,jk)=0
          end do
        end do

        do jl=1,klon
           kcbot(jl)    =  levels
           kctop(jl)    =  levels
           zqold(jl)    = 0.
           lldcum(jl)   = .false.
           resetflag(jl)= .false.
           loflag(jl)   = (.not. deepflag(jl)) .and. (levels.ge.itoppacel(jl))
        end do

! start the inner loop to search the deep convection points
      do jk=levels,2,-1
        is=0
        do jl=1,klon
         if(loflag(jl))then
            is=is+1
         endif
        enddo
        if(is.eq.0) exit

! define the variables at the departure level 
        if(jk .eq. levels) then
          do jl=1,klon
          if(loflag(jl)) then
            if((paph(jl,klev+1)-paph(jl,jk)) < 60.e2) then
              tmix=0.
              qmix=0.
              zmix=0.
              pmix=0.
              do nk=jk+2,jk,-1
              if(pmix < 50.e2) then
                dp = paph(jl,nk) - paph(jl,nk-1)
                tmix=tmix+dp*ptenh(jl,nk)
                qmix=qmix+dp*pqenh(jl,nk)
                zmix=zmix+dp*pgeoh(jl,nk)
                pmix=pmix+dp
              end if
              end do
              tmix=tmix/pmix
              qmix=qmix/pmix
              zmix=zmix/pmix
            else
              tmix=ptenh(jl,jk+1)
              qmix=pqenh(jl,jk+1)
              zmix=pgeoh(jl,jk+1)
            end if

            pqu(jl,jk+1) = qmix + deltq
            dhen(jl,jk+1)= zmix + tmix*cpd
            dh(jl,jk+1)  = dhen(jl,jk+1) + deltt*cpd
            ptu(jl,jk+1) = (dh(jl,jk+1)-pgeoh(jl,jk+1))*rcpd
            kup(jl,jk+1) = 0.5
            klab(jl,jk+1)= 1
            vptu(jl,jk+1)=ptu(jl,jk+1)*(1.+vtmpc1*pqu(jl,jk+1))
            vten(jl,jk+1)=ptenh(jl,jk+1)*(1.+vtmpc1*pqenh(jl,jk+1))
            zbuo(jl,jk+1)=(vptu(jl,jk+1)-vten(jl,jk+1))/vten(jl,jk+1)
          end if
          end do
        end if

! the next levels, we use the variables at the first level as initial values
        do jl=1,klon
           if(loflag(jl)) then
! define the fscale
             fscale = min(1.,(pqsen(jl,jk)/pqsen(jl,levels))**3)
             eta(jl) = 1.75e-3*fscale
             dz(jl)  = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
             coef(jl)= 0.5*eta(jl)*dz(jl)
             dhen(jl,jk) = pgeoh(jl,jk) + cpd*ptenh(jl,jk)
             dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk))&
     &              +(1.-coef(jl))*dh(jl,jk+1))/(1.+coef(jl))
             pqu(jl,jk) =(coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk))&
     &              +(1.-coef(jl))*pqu(jl,jk+1))/(1.+coef(jl))
             ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd
             zqold(jl) = pqu(jl,jk)
             zph(jl)=paph(jl,jk)
           end if
        end do
! check if the parcel is saturated
        ik=jk
        icall=1
        call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
 
      do jl=1,klon
        if( loflag(jl) ) then
          zdq = max((zqold(jl) - pqu(jl,jk)),0.)
          plu(jl,jk) = plu(jl,jk+1) + zdq
          zlglac=zdq*((1.-foealfa(ptu(jl,jk))) - &
                      (1.-foealfa(ptu(jl,jk+1))))
          plu(jl,jk) = 0.5*plu(jl,jk)
          dh(jl,jk) =  pgeoh(jl,jk) + cpd*(ptu(jl,jk)+ralfdcp*zlglac)
! compute the updraft speed
          vptu(jl,jk) = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))+&
                        ralfdcp*zlglac
          vten(jl,jk) = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
          zbuo(jl,jk) = (vptu(jl,jk) - vten(jl,jk))/vten(jl,jk)
          abuoy(jl,jk)=(zbuo(jl,jk)+zbuo(jl,jk+1))*0.5*g
          atop1 = 1.0 - 2.*coef(jl)
          atop2 = 2.0*dz(jl)*abuoy(jl,jk)
          abot =  1.0 + 2.*coef(jl)
          kup(jl,jk)  = (atop1*kup(jl,jk+1) + atop2) / abot
! let's find the exact cloud base
          if ( plu(jl,jk) > 0. .and. klab(jl,jk+1) == 1 ) then
              ik = jk + 1
              zqsu = foeewm(ptu(jl,ik))/paph(jl,ik)
              zqsu = min(0.5,zqsu)
              zcor = 1./(1.-vtmpc1*zqsu)
              zqsu = zqsu*zcor
              zdq = min(0.,pqu(jl,ik)-zqsu)
              zalfaw = foealfa(ptu(jl,ik))
              zfacw = c5les/((ptu(jl,ik)-c4les)**2)
              zfaci = c5ies/((ptu(jl,ik)-c4ies)**2)
              zfac = zalfaw*zfacw + (1.-zalfaw)*zfaci
              zesdp = foeewm(ptu(jl,ik))/paph(jl,ik)
              zcor = 1./(1.-vtmpc1*zesdp)
              zdqsdt = zfac*zcor*zqsu
              zdtdp = rd*ptu(jl,ik)/(cpd*paph(jl,ik))
              zdp = zdq/(zdqsdt*zdtdp)
              zcbase(jl) = paph(jl,ik) + zdp
! chose nearest half level as cloud base (jk or jk+1)
              zpdifftop = zcbase(jl) - paph(jl,jk)
              zpdiffbot = paph(jl,jk+1) - zcbase(jl)
              if ( zpdifftop > zpdiffbot .and. kup(jl,jk+1) > 0. ) then
                ikb = min(klev-1,jk+1)
                klab(jl,ikb) = 2
                klab(jl,jk) = 2
                kcbot(jl) = ikb
                plu(jl,jk+1) = 1.0e-8
              else if ( zpdifftop <= zpdiffbot .and.kup(jl,jk) > 0. ) then
                klab(jl,jk) = 2
                kcbot(jl) = jk
              end if
          end if

          if(kup(jl,jk) .lt. 0.)then
            loflag(jl) = .false.
            if(plu(jl,jk+1) .gt. 0.) then
              kctop(jl) = jk
              lldcum(jl) = .true.
            else
              lldcum(jl) = .false.
            end if
          else 
            if(plu(jl,jk) .gt. 0.)then
              klab(jl,jk)=2
            else
              klab(jl,jk)=1
            end if
          end if
        end if
      end do

      end do ! end all the levels

      needreset = .false.
      do jl=1,klon
        ikb = kcbot(jl)
        ikt = kctop(jl)
        if(paph(jl,ikb) - paph(jl,ikt) < zdnoprc) lldcum(jl) = .false.
        if(lldcum(jl)) then      
         ktype(jl)    = 1
         ldcum(jl)    = .true.
         deepflag(jl) = .true.
         wbase(jl)    = sqrt(max(2.*kup(jl,ikb),0.))
         cubot(jl)    = ikb
         cutop(jl)    = ikt
         kdpl(jl)     = levels+1
         needreset    = .true.
         resetflag(jl)= .true.
        end if
      end do

      if(needreset) then
      do jk=klev,1,-1
        do jl=1,klon
          if(resetflag(jl)) then
            ikt = kctop(jl)
            ikb = kdpl(jl)
            if(jk .le. ikb .and. jk .ge. ikt )then
             culab(jl,jk) = klab(jl,jk)
             cutu(jl,jk)  = ptu(jl,jk)
             cuqu(jl,jk)  = pqu(jl,jk)
             culu(jl,jk)  = plu(jl,jk)
            else
             culab(jl,jk) = 1
             cutu(jl,jk)  = ptenh(jl,jk)
             cuqu(jl,jk)  = pqenh(jl,jk)
             culu(jl,jk)  = 0.
            end if
            if ( jk .lt. ikt ) culab(jl,jk) = 0
          end if
        end do
      end do
      end if

      end do ! end all cycles

      return
      end subroutine cutypen

!-----------------------------------------------------------------
!    level 3 subroutines 'cuascn'
!-----------------------------------------------------------------
      subroutine cuascn &
     &    (klon,     klev,     klevp1,   klevm1,   ptenh,&
     &     pqenh,    puen,     pven,     pten,     pqen,&
     &     pqsen,    pgeo,     pgeoh,    pap,      paph,&
     &     pqte,     pverv,    klwmin,   ldcum,    phcbase,&
     &     ktype,    klab,     ptu,      pqu,      plu,&
     &     puu,      pvu,      pmfu,     pmfub,    &
     &     pmfus,    pmfuq,    pmful,    plude,    pdmfup,&
     &     kcbot,    kctop,    kctop0,   kcum,     ztmst,&
     &     pqsenh,   plglac,   lndj,     wup,      wbase,   kdpl, pmfude_rate)
      implicit none
!     this routine does the calculations for cloud ascents
!     for cumulus parameterization
!     m.tiedtke         e.c.m.w.f.     7/86 modif.  12/89
!     y.wang            iprc           11/01 modif.
!     c.zhang           iprc           05/12 modif.
!***purpose.
!   --------
!          to produce cloud ascents for cu-parametrization
!          (vertical profiles of t,q,l,u and v and corresponding
!           fluxes as well as precipitation rates)
!***interface
!   ---------
!          this routine is called from *cumastr*.
!***method.
!  --------
!          lift surface air dry-adiabatically to cloud base
!          and then calculate moist ascent for
!          entraining/detraining plume.
!          entrainment and detrainment rates differ for
!          shallow and deep cumulus convection.
!          in case there is no penetrative or shallow convection
!          check for possibility of mid level convection
!          (cloud base values calculated in *cubasmc*)
!***externals
!   ---------
!          *cuadjtqn* adjust t and q due to condensation in ascent
!          *cuentrn*  calculate entrainment/detrainment rates
!          *cubasmcn* calculate cloud base values for midlevel convection
!***reference
!   ---------
!          (tiedtke,1989)
!***input variables:
!       ptenh [ztenh] - environ temperature on half levels. (cuini)
!       pqenh [zqenh] - env. specific humidity on half levels. (cuini)
!       puen - environment wind u-component. (mssflx)
!       pven - environment wind v-component. (mssflx)
!       pten - environment temperature. (mssflx)
!       pqen - environment specific humidity. (mssflx)
!       pqsen - environment saturation specific humidity. (mssflx)
!       pgeo - geopotential. (mssflx)
!       pgeoh [zgeoh] - geopotential on half levels, (mssflx)
!       pap - pressure in pa.  (mssflx)
!       paph - pressure of half levels. (mssflx)
!       pqte - moisture convergence (delta q/delta t). (mssflx)
!       pverv - large scale vertical velocity (omega). (mssflx)
!       klwmin [ilwmin] - level of minimum omega. (cuini)
!       klab [ilab] - level label - 1: sub-cloud layer.
!                                   2: condensation level (cloud base)
!       pmfub [zmfub] - updraft mass flux at cloud base. (cumastr)
!***variables modified by cuasc:
!       ldcum - logical denoting profiles. (cubase)
!       ktype - convection type - 1: penetrative  (cumastr)
!                                 2: stratocumulus (cumastr)
!                                 3: mid-level  (cuasc)
!       ptu - cloud temperature.
!       pqu - cloud specific humidity.
!       plu - cloud liquid water (moisture condensed out)
!       puu [zuu] - cloud momentum u-component.
!       pvu [zvu] - cloud momentum v-component.
!       pmfu - updraft mass flux.
!       pmfus [zmfus] - updraft flux of dry static energy. (cubasmc)
!       pmfuq [zmfuq] - updraft flux of specific humidity.
!       pmful [zmful] - updraft flux of cloud liquid water.
!       plude - liquid water returned to environment by detrainment.
!       pdmfup [zmfup] -
!       kcbot - cloud base level. (cubase)
!       kctop - cloud top level
!       kctop0 [ictop0] - estimate of cloud top. (cumastr)
!       kcum [icum] - flag to control the call

      integer  klev,klon,klevp1,klevm1
      real(kind=kind_phys)     ptenh(klon,klev),       pqenh(klon,klev), &
     &         puen(klon,klev),        pven(klon,klev),&
     &         pten(klon,klev),        pqen(klon,klev),&
     &         pgeo(klon,klev),        pgeoh(klon,klevp1),&
     &         pap(klon,klev),         paph(klon,klevp1),&
     &         pqsen(klon,klev),       pqte(klon,klev),&
     &         pverv(klon,klev),       pqsenh(klon,klev)
      real(kind=kind_phys)     ptu(klon,klev),         pqu(klon,klev),&
     &         puu(klon,klev),         pvu(klon,klev),&
     &         pmfu(klon,klev),        zph(klon),&
     &         pmfub(klon),            &
     &         pmfus(klon,klev),       pmfuq(klon,klev),&
     &         plu(klon,klev),         plude(klon,klev),&
     &         pmful(klon,klev),       pdmfup(klon,klev)
      real(kind=kind_phys)     zdmfen(klon),           zdmfde(klon),&
     &         zmfuu(klon),            zmfuv(klon),&
     &         zpbase(klon),           zqold(klon)              
      real(kind=kind_phys)     phcbase(klon),          zluold(klon)
      real(kind=kind_phys)     zprecip(klon),          zlrain(klon,klev)
      real(kind=kind_phys)     zbuo(klon,klev),        kup(klon,klev)
      real(kind=kind_phys)     wup(klon)
      real(kind=kind_phys)     wbase(klon),            zodetr(klon,klev)
      real(kind=kind_phys)     plglac(klon,klev)

      real(kind=kind_phys)     eta(klon),dz(klon)

      integer  klwmin(klon),           ktype(klon),&
     &         klab(klon,klev),        kcbot(klon),&
     &         kctop(klon),            kctop0(klon)
      integer  lndj(klon)
      logical  ldcum(klon),            loflag(klon)
      logical  llo2,llo3,              llo1(klon)

      integer  kdpl(klon)
      real(kind=kind_phys)     zoentr(klon),           zdpmean(klon)
      real(kind=kind_phys)     pdmfen(klon,klev),      pmfude_rate(klon,klev)
! local variables
      integer  jl,jk
      integer  ikb,icum,itopm2,ik,icall,is,kcum,jlm,jll
      integer  jlx(klon)
      real(kind=kind_phys)     ztmst,zcons2,zfacbuo,zprcdgw,z_cwdrag,z_cldmax,z_cwifrac,z_cprc2
      real(kind=kind_phys)     zmftest,zmfmax,zqeen,zseen,zscde,zqude
      real(kind=kind_phys)     zmfusk,zmfuqk,zmfulk
      real(kind=kind_phys)     zbc,zbe,zkedke,zmfun,zwu,zprcon,zdt,zcbf,zzco
      real(kind=kind_phys)     zlcrit,zdfi,zc,zd,zint,zlnew,zvw,zvi,zalfaw,zrold
      real(kind=kind_phys)     zrnew,zz,zdmfeu,zdmfdu,dp
      real(kind=kind_phys)     zfac,zbuoc,zdkbuo,zdken,zvv,zarg,zchange,zxe,zxs,zdshrd
      real(kind=kind_phys)     atop1,atop2,abot
!--------------------------------
!*    1.       specify parameters
!--------------------------------
      zcons2=3./(g*ztmst)
      zfacbuo = 0.5/(1.+0.5)
      zprcdgw = cprcon*zrg
      z_cldmax = 5.e-3
      z_cwifrac = 0.5
      z_cprc2 = 0.5
      z_cwdrag = (3.0/8.0)*0.506/0.2
!---------------------------------
!     2.        set default values
!---------------------------------
      llo3 = .false.
      do jl=1,klon
        zluold(jl)=0.
        wup(jl)=0.
        zdpmean(jl)=0.
        zoentr(jl)=0.
        if(.not.ldcum(jl)) then
          ktype(jl)=0
          kcbot(jl) = -1
          pmfub(jl) = 0.
          pqu(jl,klev) = 0.
        end if
      end do

 ! initialize variout quantities     
      do jk=1,klev
      do jl=1,klon
          if(jk.ne.kcbot(jl)) plu(jl,jk)=0.
          pmfu(jl,jk)=0.
          pmfus(jl,jk)=0.
          pmfuq(jl,jk)=0.
          pmful(jl,jk)=0.
          plude(jl,jk)=0.
          plglac(jl,jk)=0.
          pdmfup(jl,jk)=0.
          zlrain(jl,jk)=0.
          zbuo(jl,jk)=0.
          kup(jl,jk)=0.
          pdmfen(jl,jk) = 0.
          pmfude_rate(jl,jk) = 0.
          if(.not.ldcum(jl).or.ktype(jl).eq.3) klab(jl,jk)=0
          if(.not.ldcum(jl).and.paph(jl,jk).lt.4.e4) kctop0(jl)=jk
      end do
      end do

      do jl = 1,klon
        if ( ktype(jl) == 3 ) ldcum(jl) = .false.
      end do
!------------------------------------------------
!     3.0      initialize values at cloud base level
!------------------------------------------------
      do jl=1,klon
        kctop(jl)=kcbot(jl)
        if(ldcum(jl)) then
          ikb = kcbot(jl)
          kup(jl,ikb) = 0.5*wbase(jl)**2
          pmfu(jl,ikb) = pmfub(jl)
          pmfus(jl,ikb) = pmfub(jl)*(cpd*ptu(jl,ikb)+pgeoh(jl,ikb))
          pmfuq(jl,ikb) = pmfub(jl)*pqu(jl,ikb)
          pmful(jl,ikb) = pmfub(jl)*plu(jl,ikb)
        end if
      end do
!
!-----------------------------------------------------------------
!     4.       do ascent: subcloud layer (klab=1) ,clouds (klab=2)
!              by doing first dry-adiabatic ascent and then
!              by adjusting t,q and l accordingly in *cuadjtqn*,
!              then check for buoyancy and set flags accordingly
!-----------------------------------------------------------------
!
      do jk=klevm1,3,-1
!             specify cloud base values for midlevel convection
!             in *cubasmc* in case there is not already convection
! ---------------------------------------------------------------------
      ik=jk
      call cubasmcn&
     &    (klon,     klev,     klevm1,   ik,      pten,&
     &     pqen,     pqsen,    puen,     pven,    pverv,&
     &     pgeo,     pgeoh,    ldcum,   ktype,   klab,  zlrain,&
     &     pmfu,     pmfub,    kcbot,   ptu,&
     &     pqu,      plu,      puu,     pvu,      pmfus,&
     &     pmfuq,    pmful,    pdmfup)
      is = 0
      jlm = 0
      do jl = 1,klon
        loflag(jl) = .false.
        zprecip(jl) = 0.
        llo1(jl) = .false.
        is = is + klab(jl,jk+1)
        if ( klab(jl,jk+1) == 0 ) klab(jl,jk) = 0
        if ( (ldcum(jl) .and. klab(jl,jk+1) == 2) .or.   &
             (ktype(jl) == 3 .and. klab(jl,jk+1) == 1) ) then
          loflag(jl) = .true.
          jlm = jlm + 1
          jlx(jlm) = jl
        end if
        zph(jl) = paph(jl,jk)
        if ( ktype(jl) == 3 .and. jk == kcbot(jl) ) then
          zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
          if ( pmfub(jl) > zmfmax ) then
            zfac = zmfmax/pmfub(jl)
            pmfu(jl,jk+1) = pmfu(jl,jk+1)*zfac
            pmfus(jl,jk+1) = pmfus(jl,jk+1)*zfac
            pmfuq(jl,jk+1) = pmfuq(jl,jk+1)*zfac
            pmfub(jl) = zmfmax
          end if
          pmfub(jl)=min(pmfub(jl),zmfmax)
        end if
      end do

      if(is.gt.0) llo3 = .true.
!
!*     specify entrainment rates in *cuentr*
! -------------------------------------
      ik=jk
      call cuentrn(klon,klev,ik,kcbot,ldcum,llo3, &
                  pgeoh,pmfu,zdmfen,zdmfde)
!
!      do adiabatic ascent for entraining/detraining plume
      if(llo3) then
! -------------------------------------------------------
!
        do jl = 1,klon
          zqold(jl) = 0.
        end do
        do jll = 1 , jlm
          jl = jlx(jll)
          zdmfde(jl) = min(zdmfde(jl),0.75*pmfu(jl,jk+1))
          if ( jk == kcbot(jl) ) then
            zoentr(jl) = -1.75e-3*(min(1.,pqen(jl,jk)/pqsen(jl,jk)) - &
                         1.)*(pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
            zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk+1)
          end if
          if ( jk < kcbot(jl) ) then
            zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
            zxs = max(pmfu(jl,jk+1)-zmfmax,0.)
            wup(jl) = wup(jl) + kup(jl,jk+1)*(pap(jl,jk+1)-pap(jl,jk))
            zdpmean(jl) = zdpmean(jl) + pap(jl,jk+1) - pap(jl,jk)
            zdmfen(jl) = zoentr(jl)
            if ( ktype(jl) >= 2 ) then
              zdmfen(jl) = 2.0*zdmfen(jl)
              zdmfde(jl) = zdmfen(jl)
            end if
            zdmfde(jl) = zdmfde(jl) * &
                         (1.6-min(1.,pqen(jl,jk)/pqsen(jl,jk)))
            zmftest = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
            zchange = max(zmftest-zmfmax,0.)
            zxe = max(zchange-zxs,0.)
            zdmfen(jl) = zdmfen(jl) - zxe
            zchange = zchange - zxe
            zdmfde(jl) = zdmfde(jl) + zchange
          end if
          pdmfen(jl,jk) = zdmfen(jl) - zdmfde(jl)
          pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
          zqeen = pqenh(jl,jk+1)*zdmfen(jl)
          zseen = (cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))*zdmfen(jl)
          zscde = (cpd*ptu(jl,jk+1)+pgeoh(jl,jk+1))*zdmfde(jl)
          zqude = pqu(jl,jk+1)*zdmfde(jl)
          plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
          zmfusk = pmfus(jl,jk+1) + zseen - zscde
          zmfuqk = pmfuq(jl,jk+1) + zqeen - zqude
          zmfulk = pmful(jl,jk+1) - plude(jl,jk)
          plu(jl,jk) = zmfulk*(1./max(cmfcmin,pmfu(jl,jk)))
          pqu(jl,jk) = zmfuqk*(1./max(cmfcmin,pmfu(jl,jk)))
          ptu(jl,jk) = (zmfusk * &
            (1./max(cmfcmin,pmfu(jl,jk)))-pgeoh(jl,jk))*rcpd
          ptu(jl,jk) = max(100.,ptu(jl,jk))
          ptu(jl,jk) = min(400.,ptu(jl,jk))
          zqold(jl) = pqu(jl,jk)
          zlrain(jl,jk) = zlrain(jl,jk+1)*(pmfu(jl,jk+1)-zdmfde(jl)) * &
                          (1./max(cmfcmin,pmfu(jl,jk)))
          zluold(jl) = plu(jl,jk)
        end do
! reset to environmental values if below departure level
        do jl = 1,klon
          if ( jk > kdpl(jl) ) then
            ptu(jl,jk) = ptenh(jl,jk)
            pqu(jl,jk) = pqenh(jl,jk)
            plu(jl,jk) = 0.
            zluold(jl) = plu(jl,jk)
          end if
        end do
!*             do corrections for moist ascent
!*             by adjusting t,q and l in *cuadjtq*
!------------------------------------------------
      ik=jk
      icall=1
!
      if ( jlm > 0 ) then
        call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
      end if
! compute the upfraft speed in cloud layer
        do jll = 1 , jlm
          jl = jlx(jll)
          if ( pqu(jl,jk) /= zqold(jl) ) then
            plglac(jl,jk) = plu(jl,jk) * &
                           ((1.-foealfa(ptu(jl,jk)))- &
                            (1.-foealfa(ptu(jl,jk+1))))
            ptu(jl,jk) = ptu(jl,jk) + ralfdcp*plglac(jl,jk)
          end if
        end do
        do jll = 1 , jlm
          jl = jlx(jll)
          if ( pqu(jl,jk) /= zqold(jl) ) then
            klab(jl,jk) = 2
            plu(jl,jk) = plu(jl,jk) + zqold(jl) - pqu(jl,jk)
            zbc = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk+1) - &
              zlrain(jl,jk+1))
            zbe = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
            zbuo(jl,jk) = zbc - zbe
! set flags for the case of midlevel convection
            if ( ktype(jl) == 3 .and. klab(jl,jk+1) == 1 ) then
              if ( zbuo(jl,jk) > -0.5 ) then
                ldcum(jl) = .true.
                kctop(jl) = jk
                kup(jl,jk) = 0.5
              else
                klab(jl,jk) = 0
                pmfu(jl,jk) = 0.
                plude(jl,jk) = 0.
                plu(jl,jk) = 0.
              end if
            end if
            if ( klab(jl,jk+1) == 2 ) then
              if ( zbuo(jl,jk) < 0. ) then
                ptenh(jl,jk) = 0.5*(pten(jl,jk)+pten(jl,jk-1))
                pqenh(jl,jk) = 0.5*(pqen(jl,jk)+pqen(jl,jk-1))
                zbuo(jl,jk) = zbc - ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
              end if
              zbuoc = (zbuo(jl,jk) / &
                (ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)))+zbuo(jl,jk+1) / &
                (ptenh(jl,jk+1)*(1.+vtmpc1*pqenh(jl,jk+1))))*0.5
              zdkbuo = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zfacbuo*zbuoc
! mixing and "pressure" gradient term in upper troposphere
              if ( zdmfen(jl) > 0. ) then
                zdken = min(1.,(1.+z_cwdrag)*zdmfen(jl) / &
                        max(cmfcmin,pmfu(jl,jk+1)))
              else
                zdken = min(1.,(1.+z_cwdrag)*zdmfde(jl) / &
                        max(cmfcmin,pmfu(jl,jk+1)))
              end if
              kup(jl,jk) = (kup(jl,jk+1)*(1.-zdken)+zdkbuo) / &
                           (1.+zdken)
              if ( zbuo(jl,jk) < 0. ) then
                zkedke = kup(jl,jk)/max(1.e-10,kup(jl,jk+1))
                zkedke = max(0.,min(1.,zkedke))
                zmfun = sqrt(zkedke)*pmfu(jl,jk+1) !* (1.6-min(1.,pqen(jl,jk) /    &
                      !  pqsen(jl,jk)))
                zdmfde(jl) = max(zdmfde(jl),pmfu(jl,jk+1)-zmfun)
                plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
                pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
              end if
              if ( zbuo(jl,jk) > -0.2  ) then
                ikb = kcbot(jl)
                zoentr(jl) = 1.75e-3*(0.3-(min(1.,pqen(jl,jk-1) /    &
                  pqsen(jl,jk-1))-1.))*(pgeoh(jl,jk-1)-pgeoh(jl,jk)) * &
                  zrg*min(1.,pqsen(jl,jk)/pqsen(jl,ikb))**3
                zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk)
              else
                zoentr(jl) = 0.
              end if
! erase values if below departure level
              if ( jk > kdpl(jl) ) then
                pmfu(jl,jk) = pmfu(jl,jk+1)
                kup(jl,jk) = 0.5
              end if
              if ( kup(jl,jk) > 0. .and. pmfu(jl,jk) > 0. ) then
                kctop(jl) = jk
                llo1(jl) = .true.
              else
                klab(jl,jk) = 0
                pmfu(jl,jk) = 0.
                kup(jl,jk) = 0.
                zdmfde(jl) = pmfu(jl,jk+1)
                plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
              end if
! save detrainment rates for updraught
              if ( pmfu(jl,jk+1) > 0. ) pmfude_rate(jl,jk) = zdmfde(jl)
            end if
          else if ( ktype(jl) == 2 .and. pqu(jl,jk) == zqold(jl) ) then
            klab(jl,jk) = 0
            pmfu(jl,jk) = 0.
            kup(jl,jk) = 0.
            zdmfde(jl) = pmfu(jl,jk+1)
            plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
            pmfude_rate(jl,jk) = zdmfde(jl)
          end if
        end do

        do jl = 1,klon
          if ( llo1(jl) ) then
! conversions only proceeds if plu is greater than a threshold liquid water
! content of 0.3 g/kg over water and 0.5 g/kg over land to prevent precipitation
! generation from small water contents.
            if ( lndj(jl).eq.1 ) then
              zdshrd = 5.e-4
            else
              zdshrd = 3.e-4
            end if
            ikb=kcbot(jl)
            if ( plu(jl,jk) > zdshrd )then
!             if ((paph(jl,ikb)-paph(jl,jk))>zdnoprc) then
              zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk+1))))
              zprcon = zprcdgw/(0.75*zwu)
! PARAMETERS FOR BERGERON-FINDEISEN PROCESS (T < -5C)
              zdt = min(rtber-rtice,max(rtber-ptu(jl,jk),0.))
              zcbf = 1. + z_cprc2*sqrt(zdt)
              zzco = zprcon*zcbf
              zlcrit = zdshrd/zcbf
              zdfi = pgeoh(jl,jk) - pgeoh(jl,jk+1)
              zc = (plu(jl,jk)-zluold(jl))
              zarg = (plu(jl,jk)/zlcrit)**2
              if ( zarg < 25.0 ) then
                zd = zzco*(1.-exp(-zarg))*zdfi
              else
                zd = zzco*zdfi
              end if
              zint = exp(-zd)
              zlnew = zluold(jl)*zint + zc/zd*(1.-zint)
              zlnew = max(0.,min(plu(jl,jk),zlnew))
              zlnew = min(z_cldmax,zlnew)
              zprecip(jl) = max(0.,zluold(jl)+zc-zlnew)
              pdmfup(jl,jk) = zprecip(jl)*pmfu(jl,jk)
              zlrain(jl,jk) = zlrain(jl,jk) + zprecip(jl)
              plu(jl,jk) = zlnew
            end if
          end if
        end do
        do jl = 1, klon
          if ( llo1(jl) ) then
            if ( zlrain(jl,jk) > 0. ) then
              zvw = 21.18*zlrain(jl,jk)**0.2
              zvi = z_cwifrac*zvw
              zalfaw = foealfa(ptu(jl,jk))
              zvv = zalfaw*zvw + (1.-zalfaw)*zvi
              zrold = zlrain(jl,jk) - zprecip(jl)
              zc = zprecip(jl)
              zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk))))
              zd = zvv/zwu
              zint = exp(-zd)
              zrnew = zrold*zint + zc/zd*(1.-zint)
              zrnew = max(0.,min(zlrain(jl,jk),zrnew))
              zlrain(jl,jk) = zrnew
            end if
          end if
        end do
        do jll = 1 , jlm
          jl = jlx(jll)
          pmful(jl,jk) = plu(jl,jk)*pmfu(jl,jk)
          pmfus(jl,jk) = (cpd*ptu(jl,jk)+pgeoh(jl,jk))*pmfu(jl,jk)
          pmfuq(jl,jk) = pqu(jl,jk)*pmfu(jl,jk)
        end do
      end if
    end do
!----------------------------------------------------------------------
! 5.       final calculations
! ------------------
      do jl = 1,klon
       if ( kctop(jl) == -1 ) ldcum(jl) = .false.
        kcbot(jl) = max(kcbot(jl),kctop(jl))
        if ( ldcum(jl) ) then
          wup(jl) = max(1.e-2,wup(jl)/max(1.,zdpmean(jl)))
          wup(jl) = sqrt(2.*wup(jl))
        end if
      end do

      return
      end subroutine cuascn
!---------------------------------------------------------
!  level 3 souroutines  
!--------------------------------------------------------
      subroutine cudlfsn   &                                                      
     &    (klon,     klev,  &                       
     &     kcbot,    kctop,    lndj,   ldcum,   &                               
     &     ptenh,    pqenh,    puen,     pven,     &                              
     &     pten,     pqsen,    pgeo,                &                             
     &     pgeoh,    paph,     ptu,      pqu,      plu,&                          
     &     puu,      pvu,      pmfub,    prfl,          &                         
     &     ptd,      pqd,      pud,      pvd,            &                        
     &     pmfd,     pmfds,    pmfdq,    pdmfdp,          &                       
     &     kdtop,    lddraf)                                                     
                                                                                 
!          this routine calculates level of free sinking for                     
!          cumulus downdrafts and specifies t,q,u and v values                   
                                                                                 
!          m.tiedtke         e.c.m.w.f.    12/86 modif. 12/89                    
                                                                                 
!          purpose.                                                              
!          --------                                                              
!          to produce lfs-values for cumulus downdrafts                          
!          for massflux cumulus parameterization                                 
                                                                                 
!          interface                                                             
!          ---------                                                             
!          this routine is called from *cumastr*.                                
!          input are environmental values of t,q,u,v,p,phi                       
!          and updraft values t,q,u and v and also                               
!          cloud base massflux and cu-precipitation rate.                        
!          it returns t,q,u and v values and massflux at lfs.                    
                                                                                 
!          method.                                                               
                                                                                 
!          check for negative buoyancy of air of equal parts of                  
!          moist environmental air and cloud air.                                
                                                                                 
!     parameter     description                                   units          
!     ---------     -----------                                   -----          
!     input parameters (integer):                                                
                                                                                 
!    *klon*         number of grid points per packet                             
!    *klev*         number of levels                                             
!    *kcbot*        cloud base level                                             
!    *kctop*        cloud top level                                              
                                                                                 
!    input parameters (logical):                                                 
                                                                                 
!    *lndj*       land sea mask (1 for land)                              
!    *ldcum*        flag: .true. for convective points                           
                                                                                 
!    input parameters (real(kind=kind_phys)):                                                    
                                                                                 
!    *ptenh*        env. temperature (t+1) on half levels          k             
!    *pqenh*        env. spec. humidity (t+1) on half levels     kg/kg           
!    *puen*         provisional environment u-velocity (t+1)      m/s            
!    *pven*         provisional environment v-velocity (t+1)      m/s            
!    *pten*         provisional environment temperature (t+1)       k            
!    *pqsen*        environment spec. saturation humidity (t+1)   kg/kg          
!    *pgeo*         geopotential                                  m2/s2          
!    *pgeoh*        geopotential on half levels                  m2/s2           
!    *paph*         provisional pressure on half levels           pa             
!    *ptu*          temperature in updrafts                        k             
!    *pqu*          spec. humidity in updrafts                   kg/kg           
!    *plu*          liquid water content in updrafts             kg/kg           
!    *puu*          u-velocity in updrafts                        m/s            
!    *pvu*          v-velocity in updrafts                        m/s            
!    *pmfub*        massflux in updrafts at cloud base           kg/(m2*s)       
                                                                                 
!    updated parameters (real(kind=kind_phys)):                                                  
                                                                                 
!    *prfl*         precipitation rate                           kg/(m2*s)       
                                                                                 
!    output parameters (real(kind=kind_phys)):                                                   
                                                                                 
!    *ptd*          temperature in downdrafts                      k             
!    *pqd*          spec. humidity in downdrafts                 kg/kg           
!    *pud*          u-velocity in downdrafts                      m/s            
!    *pvd*          v-velocity in downdrafts                      m/s            
!    *pmfd*         massflux in downdrafts                       kg/(m2*s)       
!    *pmfds*        flux of dry static energy in downdrafts       j/(m2*s)       
!    *pmfdq*        flux of spec. humidity in downdrafts         kg/(m2*s)       
!    *pdmfdp*       flux difference of precip. in downdrafts     kg/(m2*s)       
                                                                                 
!    output parameters (integer):                                                
                                                                                 
!    *kdtop*        top level of downdrafts                                      
                                                                                 
!    output parameters (logical):                                                
                                                                                 
!    *lddraf*       .true. if downdrafts exist                                   
                                                                                 
!          externals                                                             
!          ---------                                                             
!          *cuadjtq* for calculating wet bulb t and q at lfs                     
!----------------------------------------------------------------------          
      implicit none                                                       
                   
      integer  klev,klon                                                            
      real(kind=kind_phys)     ptenh(klon,klev),       pqenh(klon,klev), &                        
     &         puen(klon,klev),        pven(klon,klev),  &                        
     &         pten(klon,klev),        pqsen(klon,klev),  &                       
     &         pgeo(klon,klev),                       &                           
     &         pgeoh(klon,klev+1),     paph(klon,klev+1),&                        
     &         ptu(klon,klev),         pqu(klon,klev),   &                        
     &         puu(klon,klev),         pvu(klon,klev),   &                        
     &         plu(klon,klev),                          &                         
     &         pmfub(klon),            prfl(klon)                                
                                                                                 
      real(kind=kind_phys)     ptd(klon,klev),         pqd(klon,klev),   &                        
     &         pud(klon,klev),         pvd(klon,klev),    &                       
     &         pmfd(klon,klev),        pmfds(klon,klev),  &                       
     &         pmfdq(klon,klev),       pdmfdp(klon,klev)                         
      integer  kcbot(klon),            kctop(klon),       &                       
     &         kdtop(klon),            ikhsmin(klon)                             
      logical  ldcum(klon),                              &
     &         lddraf(klon)   
      integer  lndj(klon)                                                   
                                                                                 
      real(kind=kind_phys)     ztenwb(klon,klev),      zqenwb(klon,klev), &                       
     &         zcond(klon),            zph(klon),         &                       
     &         zhsmin(klon)                                                      
      logical  llo2(klon)                                                        
! local variables
      integer  jl,jk
      integer  is,ik,icall,ike
      real(kind=kind_phys)     zhsk,zttest,zqtest,zbuo,zmftop
                                                                                   
!----------------------------------------------------------------------          
                                                                                 
!     1.           set default values for downdrafts                             
!                  ---------------------------------                             
      do jl=1,klon                                                          
        lddraf(jl)=.false.                                                       
        kdtop(jl)=klev+1      
        ikhsmin(jl)=klev+1
        zhsmin(jl)=1.e8                                                  
      enddo                                                                      
!----------------------------------------------------------------------          
                                                                                 
!     2.           determine level of free sinking:                              
!                  downdrafts shall start at model level of minimum              
!                  of saturation moist static energy or below                    
!                  respectively                                                  
                                                                                 
!                  for every point and proceed as follows:                       
                                                                                 
!                    (1) determine level of minimum of hs                        
!                    (2) determine wet bulb environmental t and q                
!                    (3) do mixing with cumulus cloud air                        
!                    (4) check for negative buoyancy                             
!                    (5) if buoyancy>0 repeat (2) to (4) for next                
!                        level below                                             
                                                                                 
!                  the assumption is that air of downdrafts is mixture           
!                  of 50% cloud air + 50% environmental air at wet bulb          
!                  temperature (i.e. which became saturated due to               
!                  evaporation of rain and cloud water)                          
!                  ----------------------------------------------------          
      do jk=3,klev-2
         do jl=1,klon
           zhsk=cpd*pten(jl,jk)+pgeo(jl,jk) +  &
     &         foelhm(pten(jl,jk))*pqsen(jl,jk)
           if(zhsk .lt. zhsmin(jl)) then
              zhsmin(jl) = zhsk
              ikhsmin(jl)= jk
           end if
         end do
      end do


      ike=klev-3                                                                 
      do jk=3,ike                                                                
                                                                                 
!     2.1          calculate wet-bulb temperature and moisture                   
!                  for environmental air in *cuadjtq*                            
!                  -------------------------------------------                   
        is=0                                                                     
        do jl=1,klon                                                        
          ztenwb(jl,jk)=ptenh(jl,jk)                                             
          zqenwb(jl,jk)=pqenh(jl,jk)                                             
          zph(jl)=paph(jl,jk)                                                    
          llo2(jl)=ldcum(jl).and.prfl(jl).gt.0..and..not.lddraf(jl).and.   &      
     &     (jk.lt.kcbot(jl).and.jk.gt.kctop(jl)).and. jk.ge.ikhsmin(jl)                            
          if(llo2(jl))then                                                       
            is=is+1                                                              
          endif                                                                  
        enddo                                                                    
        if(is.eq.0) cycle                                                        
                                                                                 
        ik=jk                                                                    
        icall=2                                                                  
        call cuadjtqn                                           &                  
     &   ( klon, klev, ik, zph, ztenwb, zqenwb, llo2, icall)                        
                                                                                 
!     2.2          do mixing of cumulus and environmental air                    
!                  and check for negative buoyancy.                              
!                  then set values for downdraft at lfs.                         
!                  ----------------------------------------                      
        do jl=1,klon                                                        
          if(llo2(jl)) then                                                      
            zttest=0.5*(ptu(jl,jk)+ztenwb(jl,jk))                                
            zqtest=0.5*(pqu(jl,jk)+zqenwb(jl,jk))                                
            zbuo=zttest*(1.+vtmpc1  *zqtest)-                    &                  
     &       ptenh(jl,jk)*(1.+vtmpc1  *pqenh(jl,jk))                               
            zcond(jl)=pqenh(jl,jk)-zqenwb(jl,jk)                                 
            zmftop=-cmfdeps*pmfub(jl)                                            
            if(zbuo.lt.0..and.prfl(jl).gt.10.*zmftop*zcond(jl)) then             
              kdtop(jl)=jk                                                       
              lddraf(jl)=.true.                                                  
              ptd(jl,jk)=zttest                                                  
              pqd(jl,jk)=zqtest                                                  
              pmfd(jl,jk)=zmftop                                                 
              pmfds(jl,jk)=pmfd(jl,jk)*(cpd*ptd(jl,jk)+pgeoh(jl,jk))            
              pmfdq(jl,jk)=pmfd(jl,jk)*pqd(jl,jk)                                
              pdmfdp(jl,jk-1)=-0.5*pmfd(jl,jk)*zcond(jl)                         
              prfl(jl)=prfl(jl)+pdmfdp(jl,jk-1)                                  
            endif                                                                
          endif                                                                  
        enddo                                                                    
                                                                                 
      enddo                                                                      
                                                                                 
      return                                                                     
      end subroutine cudlfsn  

!---------------------------------------------------------
!  level 3 souroutines  
!--------------------------------------------------------
!**********************************************
!       subroutine cuddrafn
!**********************************************
       subroutine cuddrafn                                 &                
     &   ( klon,     klev,    lddraf                      &                                   
     &   , ptenh,    pqenh,    puen,     pven             &                
     &   , pgeo,     pgeoh,    paph,     prfl             &                
     &   , ptd,      pqd,      pud,      pvd,      pmfu   &                
     &   , pmfd,     pmfds,    pmfdq,    pdmfdp,   pmfdde_rate     )                          
                                                                          
!          this routine calculates cumulus downdraft descent              
                                                                          
!          m.tiedtke         e.c.m.w.f.    12/86 modif. 12/89             
                                                                          
!          purpose.                                                       
!          --------                                                       
!          to produce the vertical profiles for cumulus downdrafts        
!          (i.e. t,q,u and v and fluxes)                                  
                                                                          
!          interface                                                      
!          ---------                                                      
                                                                          
!          this routine is called from *cumastr*.                         
!          input is t,q,p,phi,u,v at half levels.                         
!          it returns fluxes of s,q and evaporation rate                  
!          and u,v at levels where downdraft occurs                       
                                                                          
!          method.                                                        
!          --------                                                       
!          calculate moist descent for entraining/detraining plume by     
!          a) moving air dry-adiabatically to next level below and        
!          b) correcting for evaporation to obtain saturated state.       
                                                                          
!     parameter     description                                   units   
!     ---------     -----------                                   -----   
!     input parameters (integer):                                         
                                                                          
!    *klon*         number of grid points per packet                      
!    *klev*         number of levels                                      
                                                                          
!    input parameters (logical):                                          
                                                                          
!    *lddraf*       .true. if downdrafts exist                            
                                                                          
!    input parameters (real(kind=kind_phys)):                                             
                                                                          
!    *ptenh*        env. temperature (t+1) on half levels          k      
!    *pqenh*        env. spec. humidity (t+1) on half levels     kg/kg    
!    *puen*         provisional environment u-velocity (t+1)      m/s     
!    *pven*         provisional environment v-velocity (t+1)      m/s     
!    *pgeo*         geopotential                                  m2/s2   
!    *pgeoh*        geopotential on half levels                  m2/s2    
!    *paph*         provisional pressure on half levels           pa      
!    *pmfu*         massflux updrafts                           kg/(m2*s) 
                                                                          
!    updated parameters (real(kind=kind_phys)):                                           
                                                                          
!    *prfl*         precipitation rate                           kg/(m2*s)
                                                                          
!    output parameters (real(kind=kind_phys)):                                            
                                                                          
!    *ptd*          temperature in downdrafts                      k      
!    *pqd*          spec. humidity in downdrafts                 kg/kg    
!    *pud*          u-velocity in downdrafts                      m/s     
!    *pvd*          v-velocity in downdrafts                      m/s     
!    *pmfd*         massflux in downdrafts                       kg/(m2*s)
!    *pmfds*        flux of dry static energy in downdrafts       j/(m2*s)
!    *pmfdq*        flux of spec. humidity in downdrafts         kg/(m2*s)
!    *pdmfdp*       flux difference of precip. in downdrafts     kg/(m2*s)
                                                                          
!          externals                                                      
!          ---------                                                      
!          *cuadjtq* for adjusting t and q due to evaporation in          
!          saturated descent                                              
!----------------------------------------------------------------------   
      implicit none
      
      integer  klev,klon                                                             
      real(kind=kind_phys)     ptenh(klon,klev),       pqenh(klon,klev),   &               
     &         puen(klon,klev),        pven(klon,klev),    &               
     &         pgeoh(klon,klev+1),     paph(klon,klev+1),  &               
     &         pgeo(klon,klev),        pmfu(klon,klev)                    
                                                                          
      real(kind=kind_phys)     ptd(klon,klev),         pqd(klon,klev),     &               
     &         pud(klon,klev),         pvd(klon,klev),     &               
     &         pmfd(klon,klev),        pmfds(klon,klev),   &               
     &         pmfdq(klon,klev),       pdmfdp(klon,klev),  &               
     &         prfl(klon)     
      real(kind=kind_phys)     pmfdde_rate(klon,klev)                                            
      logical  lddraf(klon)                                               
                                                                          
      real(kind=kind_phys)     zdmfen(klon),           zdmfde(klon),       &               
     &         zcond(klon),            zoentr(klon),       &               
     &         zbuoy(klon)                                                
      real(kind=kind_phys)     zph(klon)                         
      logical  llo2(klon)                                   
      logical  llo1
! local variables
      integer  jl,jk
      integer  is,ik,icall,ike, itopde(klon)
      real(kind=kind_phys)     zentr,zdz,zzentr,zseen,zqeen,zsdde,zqdde,zdmfdp
      real(kind=kind_phys)     zmfdsk,zmfdqk,zbuo,zrain,zbuoyz,zmfduk,zmfdvk
                                                                
!----------------------------------------------------------------------   
!     1.           calculate moist descent for cumulus downdraft by       
!                     (a) calculating entrainment/detrainment rates,      
!                         including organized entrainment dependent on    
!                         negative buoyancy and assuming                  
!                         linear decrease of massflux in pbl              
!                     (b) doing moist descent - evaporative cooling       
!                         and moistening is calculated in *cuadjtq*       
!                     (c) checking for negative buoyancy and              
!                         specifying final t,q,u,v and downward fluxes    
!                    -------------------------------------------------    
      do jl=1,klon                                                   
        zoentr(jl)=0.                                                     
        zbuoy(jl)=0.                                                      
        zdmfen(jl)=0.                                                     
        zdmfde(jl)=0.
      enddo 

      do jk=klev,1,-1
       do jl=1,klon
         pmfdde_rate(jl,jk) = 0.
         if((paph(jl,klev+1)-paph(jl,jk)).lt. 60.e2) itopde(jl) = jk
       end do
      end do                                                              
                                                                 
      do jk=3,klev                                                        
        is=0                                                              
        do jl=1,klon                                                 
          zph(jl)=paph(jl,jk)                                             
          llo2(jl)=lddraf(jl).and.pmfd(jl,jk-1).lt.0.                     
          if(llo2(jl)) then                                               
            is=is+1                                                       
          endif                                                           
        enddo                                                             
        if(is.eq.0) cycle                                                 
                                                                          
        do jl=1,klon                                                 
          if(llo2(jl)) then    
            zentr = entrdd*pmfd(jl,jk-1)*(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg
            zdmfen(jl)=zentr                                              
            zdmfde(jl)=zentr                                              
          endif                                                           
        enddo 
                                                            
        do jl=1,klon
          if(llo2(jl)) then
          if(jk.gt.itopde(jl)) then
            zdmfen(jl)=0.
            zdmfde(jl)=pmfd(jl,itopde(jl))*                    &
     &       (paph(jl,jk)-paph(jl,jk-1))/                  &
     &       (paph(jl,klev+1)-paph(jl,itopde(jl)))
          endif
          endif
        enddo

        do jl=1,klon
          if(llo2(jl)) then
          if(jk.le.itopde(jl)) then
            zdz=-(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg
            zzentr=zoentr(jl)*zdz*pmfd(jl,jk-1)
            zdmfen(jl)=zdmfen(jl)+zzentr
            zdmfen(jl)=max(zdmfen(jl),0.3*pmfd(jl,jk-1))
            zdmfen(jl)=max(zdmfen(jl),-0.75*pmfu(jl,jk)-   &
   &         (pmfd(jl,jk-1)-zdmfde(jl)))
            zdmfen(jl)=min(zdmfen(jl),0.)
          endif
          endif
        enddo

        do jl=1,klon                                                 
          if(llo2(jl)) then                                               
            pmfd(jl,jk)=pmfd(jl,jk-1)+zdmfen(jl)-zdmfde(jl)
            zseen=(cpd*ptenh(jl,jk-1)+pgeoh(jl,jk-1))*zdmfen(jl)         
            zqeen=pqenh(jl,jk-1)*zdmfen(jl)                               
            zsdde=(cpd*ptd(jl,jk-1)+pgeoh(jl,jk-1))*zdmfde(jl)           
            zqdde=pqd(jl,jk-1)*zdmfde(jl)                                 
            zmfdsk=pmfds(jl,jk-1)+zseen-zsdde                             
            zmfdqk=pmfdq(jl,jk-1)+zqeen-zqdde                             
            pqd(jl,jk)=zmfdqk*(1./min(-cmfcmin,pmfd(jl,jk)))              
            ptd(jl,jk)=(zmfdsk*(1./min(-cmfcmin,pmfd(jl,jk)))-&            
     &                  pgeoh(jl,jk))*rcpd                                
            ptd(jl,jk)=min(400.,ptd(jl,jk))                               
            ptd(jl,jk)=max(100.,ptd(jl,jk))                               
            zcond(jl)=pqd(jl,jk)                                          
          endif                                                           
        enddo                                                             
                                                                          
        ik=jk                                                             
        icall=2                                                           
        call cuadjtqn(klon, klev, ik, zph, ptd, pqd, llo2, icall )                
                                                                          
        do jl=1,klon                                                 
          if(llo2(jl)) then                                               
            zcond(jl)=zcond(jl)-pqd(jl,jk)
            zbuo=ptd(jl,jk)*(1.+vtmpc1  *pqd(jl,jk))-          &             
     &      ptenh(jl,jk)*(1.+vtmpc1  *pqenh(jl,jk))                         
            if(prfl(jl).gt.0..and.pmfu(jl,jk).gt.0.) then                 
              zrain=prfl(jl)/pmfu(jl,jk)                                  
              zbuo=zbuo-ptd(jl,jk)*zrain
            endif                                                         
            if(zbuo.ge.0 .or. prfl(jl).le.(pmfd(jl,jk)*zcond(jl))) then    
              pmfd(jl,jk)=0.
              zbuo=0.                                              
            endif
            pmfds(jl,jk)=(cpd*ptd(jl,jk)+pgeoh(jl,jk))*pmfd(jl,jk)       
            pmfdq(jl,jk)=pqd(jl,jk)*pmfd(jl,jk)                           
            zdmfdp=-pmfd(jl,jk)*zcond(jl)                                 
            pdmfdp(jl,jk-1)=zdmfdp                                        
            prfl(jl)=prfl(jl)+zdmfdp                                      
                                                                          
! compute organized entrainment for use at next level                     
            zbuoyz=zbuo/ptenh(jl,jk)                                      
            zbuoyz=min(zbuoyz,0.0)                                        
            zdz=-(pgeo(jl,jk-1)-pgeo(jl,jk))
            zbuoy(jl)=zbuoy(jl)+zbuoyz*zdz
            zoentr(jl)=g*zbuoyz*0.5/(1.+zbuoy(jl)) 
            pmfdde_rate(jl,jk) = -zdmfde(jl)
          endif                                                           
        enddo
                                                             
      enddo                                                               
                                                                          
      return                                                              
      end subroutine cuddrafn
!---------------------------------------------------------
!  level 3 souroutines  
!--------------------------------------------------------
       subroutine cuflxn         &                                                 
     &  (  klon,     klev,     ztmst &                                                             
     &  ,  pten,     pqen,     pqsen,    ptenh,    pqenh   &                    
     &  ,  paph,     pap,      pgeoh,    lndj,   ldcum   &                    
     &  ,  kcbot,    kctop,    kdtop,    ktopm2            &                    
     &  ,  ktype,    lddraf                                &                    
     &  ,  pmfu,     pmfd,     pmfus,    pmfds             &                    
     &  ,  pmfuq,    pmfdq,    pmful,    plude             &                    
     &  ,  pdmfup,   pdmfdp,   pdpmel,   plglac            &                    
     &  ,  prain,    pmfdde_rate, pmflxr, pmflxs )                                         
                                                                               
!          m.tiedtke         e.c.m.w.f.     7/86 modif. 12/89                  
                                                                               
!          purpose                                                             
!          -------                                                             
                                                                               
!          this routine does the final calculation of convective               
!          fluxes in the cloud layer and in the subcloud layer                 
                                                                               
!          interface                                                           
!          ---------                                                           
!          this routine is called from *cumastr*.                              
                                                                               
                                                                               
!     parameter     description                                   units        
!     ---------     -----------                                   -----        
!     input parameters (integer):                                              
                                                                               
!    *klon*         number of grid points per packet                           
!    *klev*         number of levels                                           
!    *kcbot*        cloud base level                                           
!    *kctop*        cloud top level                                            
!    *kdtop*        top level of downdrafts                                    
                                                                               
!    input parameters (logical):                                               
                                                                               
!    *lndj*       land sea mask (1 for land)                            
!    *ldcum*        flag: .true. for convective points                         
                                                                               
!    input parameters (real(kind=kind_phys)):                                                  
                                                                               
!    *ztmst*        time step for the physics                       s          
!    *pten*         provisional environment temperature (t+1)       k          
!    *pqen*         provisional environment spec. humidity (t+1)  kg/kg        
!    *pqsen*        environment spec. saturation humidity (t+1)   kg/kg        
!    *ptenh*        env. temperature (t+1) on half levels           k          
!    *pqenh*        env. spec. humidity (t+1) on half levels      kg/kg        
!    *paph*         provisional pressure on half levels            pa          
!    *pap*          provisional pressure on full levels            pa          
!    *pgeoh*        geopotential on half levels                   m2/s2        
                                                                               
!    updated parameters (integer):                                             
                                                                               
!    *ktype*        set to zero if ldcum=.false.                               
                                                                               
!    updated parameters (logical):                                             
                                                                               
!    *lddraf*       set to .false. if ldcum=.false. or kdtop<kctop             
                                                                               
!    updated parameters (real(kind=kind_phys)):                                                
                                                                               
!    *pmfu*         massflux in updrafts                          kg/(m2*s)    
!    *pmfd*         massflux in downdrafts                        kg/(m2*s)    
!    *pmfus*        flux of dry static energy in updrafts          j/(m2*s)    
!    *pmfds*        flux of dry static energy in downdrafts        j/(m2*s)    
!    *pmfuq*        flux of spec. humidity in updrafts            kg/(m2*s)    
!    *pmfdq*        flux of spec. humidity in downdrafts          kg/(m2*s)    
!    *pmful*        flux of liquid water in updrafts              kg/(m2*s)    
!    *plude*        detrained liquid water                        kg/(m3*s)    
!    *pdmfup*       flux difference of precip. in updrafts        kg/(m2*s)    
!    *pdmfdp*       flux difference of precip. in downdrafts      kg/(m2*s)    
                                                                               
!    output parameters (real(kind=kind_phys)):                                                 
                                                                               
!    *pdpmel*       change in precip.-fluxes due to melting       kg/(m2*s)    
!    *plglac*       flux of frozen cloud water in updrafts        kg/(m2*s)    
!    *pmflxr*       convective rain flux                          kg/(m2*s)    
!    *pmflxs*       convective snow flux                          kg/(m2*s)    
!    *prain*        total precip. produced in conv. updrafts      kg/(m2*s)    
!                   (no evaporation in downdrafts)                             
                                                                               
!          externals                                                           
!          ---------                                                           
!          none                                                                
!----------------------------------------------------------------------        
      implicit none                                                     

      integer  klev,klon,ktopm2                                                                         
      real(kind=kind_phys)     pten(klon,klev),        ptenh(klon,klev),          &             
     &         pqen(klon,klev),        pqsen(klon,klev),          &             
     &         pqenh(klon,klev),       pap(klon,klev),            &             
     &         paph(klon,klev+1),      pgeoh(klon,klev+1),        &             
     &         plglac(klon,klev)                                               
                                                                               
      real(kind=kind_phys)     pmfu(klon,klev),        pmfd(klon,klev),           &             
     &         pmfus(klon,klev),       pmfds(klon,klev),          &             
     &         pmfuq(klon,klev),       pmfdq(klon,klev),          &             
     &         pdmfup(klon,klev),      pdmfdp(klon,klev),         &             
     &         pdpmel(klon,klev),      prain(klon),               &             
     &         pmful(klon,klev),       plude(klon,klev),          &             
     &         pmflxr(klon,klev+1),    pmflxs(klon,klev+1)  
      real(kind=kind_phys)     pmfdde_rate(klon,klev)               
      integer  kcbot(klon),            kctop(klon),               &             
     &         kdtop(klon),            ktype(klon)                             
      logical  lddraf(klon),                                &
     &         ldcum(klon)  
      integer  lndj(klon)
! local variables
      integer  jl,jk
      integer  is,ik,icall,ike,ikb
      real(kind=kind_phys)     ztmst,ztaumel,zcons1a,zcons1,zcons2,zcucov,zcpecons
      real(kind=kind_phys)     zalfaw,zrfl,zdrfl1,zrnew,zrmin,zrfln,zdrfl,zdenom
      real(kind=kind_phys)     zpdr,zpds,zzp,zfac,zsnmlt
      real(kind=kind_phys)     rhevap(klon)
      integer  idbas(klon)
      logical  llddraf
!--------------------------------------------------------------------          
!*             specify constants  

      ztaumel=18000.
      zcons1a=cpd/(alf*g*ztaumel)
      zcons2=3./(g*ztmst)
      zcucov=0.05
      zcpecons=5.44e-4/g

!*    1.0          determine final convective fluxes                           
!                  ---------------------------------                           
      do jl=1,klon
        prain(jl)=0.
        if(.not.ldcum(jl).or.kdtop(jl).lt.kctop(jl)) lddraf(jl)=.false.
        if(.not.ldcum(jl)) ktype(jl)=0
        idbas(jl) = klev
        if(lndj(jl) .eq. 1) then
          rhevap(jl) = 0.7
        else
          rhevap(jl) = 0.9
        end if
      enddo

      do jk=ktopm2,klev
        ikb = min(jk+1,klev)
        do jl=1,klon
          pmflxr(jl,jk) = 0.
          pmflxs(jl,jk) = 0.
          pdpmel(jl,jk) = 0.
          if(ldcum(jl).and.jk.ge.kctop(jl)) then
            pmfus(jl,jk)=pmfus(jl,jk)-pmfu(jl,jk)*           &
     &       (cpd*ptenh(jl,jk)+pgeoh(jl,jk))
            pmfuq(jl,jk)=pmfuq(jl,jk)-pmfu(jl,jk)*pqenh(jl,jk)
            plglac(jl,jk)=pmfu(jl,jk)*plglac(jl,jk)
            llddraf = lddraf(jl) .and. jk >= kdtop(jl)
            if ( llddraf .and.jk.ge.kdtop(jl)) then
              pmfds(jl,jk) = pmfds(jl,jk)-pmfd(jl,jk) * &
                           (cpd*ptenh(jl,jk)+pgeoh(jl,jk))
              pmfdq(jl,jk) = pmfdq(jl,jk)-pmfd(jl,jk)*pqenh(jl,jk)
            else
              pmfd(jl,jk) = 0.
              pmfds(jl,jk) = 0.
              pmfdq(jl,jk) = 0.
              pdmfdp(jl,jk-1) = 0.
            end if
            if ( llddraf .and. pmfd(jl,jk) < 0. .and. &
               abs(pmfd(jl,ikb)) < 1.e-20 ) then
               idbas(jl) = jk
            end if
          else
            pmfu(jl,jk)=0.
            pmfd(jl,jk)=0.
            pmfus(jl,jk)=0.
            pmfds(jl,jk)=0.
            pmfuq(jl,jk)=0.
            pmfdq(jl,jk)=0.
            pmful(jl,jk)=0.
            plglac(jl,jk)=0.
            pdmfup(jl,jk-1)=0.
            pdmfdp(jl,jk-1)=0.
            plude(jl,jk-1)=0.
          endif
        enddo
      enddo

      do jl=1,klon
        pmflxr(jl,klev+1) = 0.
        pmflxs(jl,klev+1) = 0.
      end do
      do jl=1,klon
        if(ldcum(jl)) then
          ikb=kcbot(jl)
          ik=ikb+1
          zzp=((paph(jl,klev+1)-paph(jl,ik))/                  &
     &     (paph(jl,klev+1)-paph(jl,ikb)))
          if(ktype(jl).eq.3) then
            zzp=zzp**2
          endif
          pmfu(jl,ik)=pmfu(jl,ikb)*zzp
          pmfus(jl,ik)=(pmfus(jl,ikb)-                         &
     &         foelhm(ptenh(jl,ikb))*pmful(jl,ikb))*zzp
          pmfuq(jl,ik)=(pmfuq(jl,ikb)+pmful(jl,ikb))*zzp
          pmful(jl,ik)=0.
        endif
      enddo

      do jk=ktopm2,klev
        do jl=1,klon
          if(ldcum(jl).and.jk.gt.kcbot(jl)+1) then
            ikb=kcbot(jl)+1
            zzp=((paph(jl,klev+1)-paph(jl,jk))/                &
     &       (paph(jl,klev+1)-paph(jl,ikb)))
            if(ktype(jl).eq.3) then
              zzp=zzp**2
            endif
            pmfu(jl,jk)=pmfu(jl,ikb)*zzp
            pmfus(jl,jk)=pmfus(jl,ikb)*zzp
            pmfuq(jl,jk)=pmfuq(jl,ikb)*zzp
            pmful(jl,jk)=0.
          endif
          ik = idbas(jl)
          llddraf = lddraf(jl) .and. jk > ik .and. ik < klev
          if ( llddraf .and. ik == kcbot(jl)+1 ) then
           zzp = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ik)))
           if ( ktype(jl) == 3 ) zzp = zzp*zzp
            pmfd(jl,jk) = pmfd(jl,ik)*zzp
            pmfds(jl,jk) = pmfds(jl,ik)*zzp
            pmfdq(jl,jk) = pmfdq(jl,ik)*zzp
            pmfdde_rate(jl,jk) = -(pmfd(jl,jk-1)-pmfd(jl,jk))
           else if ( llddraf .and. ik /= kcbot(jl)+1 .and. jk == ik+1 ) then
            pmfdde_rate(jl,jk) = -(pmfd(jl,jk-1)-pmfd(jl,jk))
           end if
        enddo
      enddo
!*    2.            calculate rain/snow fall rates                             
!*                  calculate melting of snow                                  
!*                  calculate evaporation of precip                            
!                   -------------------------------                            

        do jk=ktopm2,klev
        do jl=1,klon
          if(ldcum(jl) .and. jk >=kctop(jl)-1 ) then
            prain(jl)=prain(jl)+pdmfup(jl,jk)
            if(pmflxs(jl,jk).gt.0..and.pten(jl,jk).gt.tmelt) then
              zcons1=zcons1a*(1.+0.5*(pten(jl,jk)-tmelt))
              zfac=zcons1*(paph(jl,jk+1)-paph(jl,jk))
              zsnmlt=min(pmflxs(jl,jk),zfac*(pten(jl,jk)-tmelt))
              pdpmel(jl,jk)=zsnmlt
              pqsen(jl,jk)=foeewm(pten(jl,jk)-zsnmlt/zfac)/pap(jl,jk)
            endif
            zalfaw=foealfa(pten(jl,jk))
          !
          ! No liquid precipitation above melting level
          !
            if ( pten(jl,jk) < tmelt .and. zalfaw > 0. ) then
              plglac(jl,jk) = plglac(jl,jk)+zalfaw*(pdmfup(jl,jk)+pdmfdp(jl,jk))
              zalfaw = 0.
            end if
            pmflxr(jl,jk+1)=pmflxr(jl,jk)+zalfaw*               &
     &       (pdmfup(jl,jk)+pdmfdp(jl,jk))+pdpmel(jl,jk)
            pmflxs(jl,jk+1)=pmflxs(jl,jk)+(1.-zalfaw)*          &
     &       (pdmfup(jl,jk)+pdmfdp(jl,jk))-pdpmel(jl,jk)
            if(pmflxr(jl,jk+1)+pmflxs(jl,jk+1).lt.0.0) then
              pdmfdp(jl,jk)=-(pmflxr(jl,jk)+pmflxs(jl,jk)+pdmfup(jl,jk))
              pmflxr(jl,jk+1)=0.0
              pmflxs(jl,jk+1)=0.0
              pdpmel(jl,jk)  =0.0
            else if ( pmflxr(jl,jk+1) < 0. ) then
              pmflxs(jl,jk+1) = pmflxs(jl,jk+1)+pmflxr(jl,jk+1)
              pmflxr(jl,jk+1) = 0.
            else if ( pmflxs(jl,jk+1) < 0. ) then
              pmflxr(jl,jk+1) = pmflxr(jl,jk+1)+pmflxs(jl,jk+1)
              pmflxs(jl,jk+1) = 0.
            end if
          endif
        enddo
      enddo
      do jk=ktopm2,klev
        do jl=1,klon
          if(ldcum(jl).and.jk.ge.kcbot(jl)) then
            zrfl=pmflxr(jl,jk)+pmflxs(jl,jk)
            if(zrfl.gt.1.e-20) then
              zdrfl1=zcpecons*max(0.,pqsen(jl,jk)-pqen(jl,jk))*zcucov* &
     &         (sqrt(paph(jl,jk)/paph(jl,klev+1))/5.09e-3*             &
     &         zrfl/zcucov)**0.5777*                                   &
     &         (paph(jl,jk+1)-paph(jl,jk))
              zrnew=zrfl-zdrfl1
              zrmin=zrfl-zcucov*max(0.,rhevap(jl)*pqsen(jl,jk) &
     &              -pqen(jl,jk)) *zcons2*(paph(jl,jk+1)-paph(jl,jk))
              zrnew=max(zrnew,zrmin)
              zrfln=max(zrnew,0.)
              zdrfl=min(0.,zrfln-zrfl)
              zdenom=1./max(1.e-20,pmflxr(jl,jk)+pmflxs(jl,jk))
              zalfaw=foealfa(pten(jl,jk))
              if ( pten(jl,jk) < tmelt ) zalfaw = 0.
              zpdr=zalfaw*pdmfdp(jl,jk)
              zpds=(1.-zalfaw)*pdmfdp(jl,jk)
              pmflxr(jl,jk+1)=pmflxr(jl,jk)+zpdr         &
     &         +pdpmel(jl,jk)+zdrfl*pmflxr(jl,jk)*zdenom
              pmflxs(jl,jk+1)=pmflxs(jl,jk)+zpds         &
     &         -pdpmel(jl,jk)+zdrfl*pmflxs(jl,jk)*zdenom
              pdmfup(jl,jk)=pdmfup(jl,jk)+zdrfl
              if ( pmflxr(jl,jk+1)+pmflxs(jl,jk+1) < 0. ) then
                pdmfup(jl,jk) = pdmfup(jl,jk)-(pmflxr(jl,jk+1)+pmflxs(jl,jk+1))
                pmflxr(jl,jk+1) = 0.
                pmflxs(jl,jk+1) = 0.
                pdpmel(jl,jk)   = 0.
              else if ( pmflxr(jl,jk+1) < 0. ) then
                pmflxs(jl,jk+1) = pmflxs(jl,jk+1)+pmflxr(jl,jk+1)
                pmflxr(jl,jk+1) = 0.
              else if ( pmflxs(jl,jk+1) < 0. ) then
                pmflxr(jl,jk+1) = pmflxr(jl,jk+1)+pmflxs(jl,jk+1)
                pmflxs(jl,jk+1) = 0.
              end if
            else
              pmflxr(jl,jk+1)=0.0
              pmflxs(jl,jk+1)=0.0
              pdmfdp(jl,jk)=0.0
              pdpmel(jl,jk)=0.0
            endif
          endif
        enddo
      enddo
                                      
      return
      end subroutine cuflxn 
!---------------------------------------------------------
!  level 3 souroutines  
!--------------------------------------------------------
     subroutine cudtdqn(klon,klev,ktopm2,kctop,kdtop,ldcum, &
                     lddraf,ztmst,paph,pgeoh,pgeo,pten,ptenh,pqen,  &
                     pqenh,pqsen,plglac,plude,pmfu,pmfd,pmfus,pmfds, &
                     pmfuq,pmfdq,pmful,pdmfup,pdmfdp,pdpmel,ptent,ptenq,pcte)
    implicit none
    integer  klon,klev,ktopm2
    integer  kctop(klon),  kdtop(klon)
    logical  ldcum(klon),  lddraf(klon)
    real(kind=kind_phys)     ztmst
    real(kind=kind_phys)     paph(klon,klev+1), pgeoh(klon,klev+1)
    real(kind=kind_phys)     pgeo(klon,klev),   pten(klon,klev), &
             pqen(klon,klev),   ptenh(klon,klev),&
             pqenh(klon,klev),  pqsen(klon,klev),&
             plglac(klon,klev), plude(klon,klev)
    real(kind=kind_phys)     pmfu(klon,klev),   pmfd(klon,klev),&
             pmfus(klon,klev),  pmfds(klon,klev),&
             pmfuq(klon,klev),  pmfdq(klon,klev),&
             pmful(klon,klev),  pdmfup(klon,klev),&
             pdpmel(klon,klev), pdmfdp(klon,klev)
    real(kind=kind_phys)     ptent(klon,klev),  ptenq(klon,klev)
    real(kind=kind_phys)     pcte(klon,klev)

! local variables
    integer  jk , ik , jl
    real(kind=kind_phys)     zalv , zzp
    real(kind=kind_phys)     zdtdt(klon,klev) , zdqdt(klon,klev) , zdp(klon,klev)
    !*    1.0          SETUP AND INITIALIZATIONS
    ! -------------------------
    do jk = 1 , klev
      do jl = 1, klon
        if ( ldcum(jl) ) then
          zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk))
        end if
      end do
    end do

    !-----------------------------------------------------------------------
    !*    2.0          COMPUTE TENDENCIES
    ! ------------------
    do jk = ktopm2 , klev
      if ( jk < klev ) then
        do jl = 1,klon
          if ( ldcum(jl) ) then
            zalv = foelhm(pten(jl,jk))
            zdtdt(jl,jk) = zdp(jl,jk)*rcpd * &
              (pmfus(jl,jk+1)-pmfus(jl,jk)+pmfds(jl,jk+1) - &
               pmfds(jl,jk)+alf*plglac(jl,jk)-alf*pdpmel(jl,jk) - &
               zalv*(pmful(jl,jk+1)-pmful(jl,jk)-plude(jl,jk)-pdmfup(jl,jk)-pdmfdp(jl,jk)))
            zdqdt(jl,jk) = zdp(jl,jk)*(pmfuq(jl,jk+1) - &
              pmfuq(jl,jk)+pmfdq(jl,jk+1)-pmfdq(jl,jk)+pmful(jl,jk+1) - &
              pmful(jl,jk)-plude(jl,jk)-pdmfup(jl,jk)-pdmfdp(jl,jk))
          end if
        end do
      else
        do jl = 1,klon
          if ( ldcum(jl) ) then
            zalv = foelhm(pten(jl,jk))
            zdtdt(jl,jk) = -zdp(jl,jk)*rcpd * &
              (pmfus(jl,jk)+pmfds(jl,jk)+alf*pdpmel(jl,jk) - &
               zalv*(pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)+plude(jl,jk)))
            zdqdt(jl,jk) = -zdp(jl,jk)*(pmfuq(jl,jk) + plude(jl,jk) + &
              pmfdq(jl,jk)+(pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)))
          end if
        end do
      end if
    end do
  !---------------------------------------------------------------
  !*  3.0          UPDATE TENDENCIES
  !   -----------------
    do jk = ktopm2 , klev
     do jl = 1, klon
       if ( ldcum(jl) ) then
         ptent(jl,jk) = ptent(jl,jk) + zdtdt(jl,jk)
         ptenq(jl,jk) = ptenq(jl,jk) + zdqdt(jl,jk)
         pcte(jl,jk)  = zdp(jl,jk)*plude(jl,jk)
       end if
     end do
    end do

    return
  end subroutine cudtdqn
!---------------------------------------------------------
!  level 3 souroutines  
!--------------------------------------------------------
    subroutine cududvn(klon,klev,ktopm2,ktype,kcbot,kctop,ldcum,  &
                    ztmst,paph,puen,pven,pmfu,pmfd,puu,pud,pvu,pvd,ptenu, &
                    ptenv)
    implicit none
    integer  klon,klev,ktopm2
    integer  ktype(klon), kcbot(klon), kctop(klon)
    logical  ldcum(klon)
    real(kind=kind_phys)     ztmst
    real(kind=kind_phys)     paph(klon,klev+1)
    real(kind=kind_phys)     puen(klon,klev),  pven(klon,klev),&
             pmfu(klon,klev),  pmfd(klon,klev),&
             puu(klon,klev),   pud(klon,klev),&
             pvu(klon,klev),   pvd(klon,klev)
    real(kind=kind_phys)     ptenu(klon,klev), ptenv(klon,klev)

!local variables
    real(kind=kind_phys)     zuen(klon,klev) , zven(klon,klev) , zmfuu(klon,klev), &
             zmfdu(klon,klev), zmfuv(klon,klev), zmfdv(klon,klev)

    integer  ik , ikb , jk , jl
    real(kind=kind_phys)     zzp, zdtdt
   
    real(kind=kind_phys)     zdudt(klon,klev), zdvdt(klon,klev), zdp(klon,klev)
!
    do jk = 1 , klev
      do jl = 1, klon
        if ( ldcum(jl) ) then
          zuen(jl,jk) = puen(jl,jk)
          zven(jl,jk) = pven(jl,jk)
          zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk))
        end if
      end do
    end do
!----------------------------------------------------------------------
!*    1.0          CALCULATE FLUXES AND UPDATE U AND V TENDENCIES
! ----------------------------------------------
    do jk = ktopm2 , klev
      ik = jk - 1
      do jl = 1,klon
        if ( ldcum(jl) ) then
          zmfuu(jl,jk) = pmfu(jl,jk)*(puu(jl,jk)-zuen(jl,ik))
          zmfuv(jl,jk) = pmfu(jl,jk)*(pvu(jl,jk)-zven(jl,ik))
          zmfdu(jl,jk) = pmfd(jl,jk)*(pud(jl,jk)-zuen(jl,ik))
          zmfdv(jl,jk) = pmfd(jl,jk)*(pvd(jl,jk)-zven(jl,ik))
        end if
      end do
    end do
    ! linear fluxes below cloud
      do jk = ktopm2 , klev
        do jl = 1, klon
          if ( ldcum(jl) .and. jk > kcbot(jl) ) then
            ikb = kcbot(jl)
            zzp = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ikb)))
            if ( ktype(jl) == 3 ) zzp = zzp*zzp
            zmfuu(jl,jk) = zmfuu(jl,ikb)*zzp
            zmfuv(jl,jk) = zmfuv(jl,ikb)*zzp
            zmfdu(jl,jk) = zmfdu(jl,ikb)*zzp
            zmfdv(jl,jk) = zmfdv(jl,ikb)*zzp
          end if
        end do
      end do
!----------------------------------------------------------------------
!*    2.0          COMPUTE TENDENCIES
! ------------------
    do jk = ktopm2 , klev
      if ( jk < klev ) then
        ik = jk + 1
        do jl = 1,klon
          if ( ldcum(jl) ) then
            zdudt(jl,jk) = zdp(jl,jk) * &
                          (zmfuu(jl,ik)-zmfuu(jl,jk)+zmfdu(jl,ik)-zmfdu(jl,jk))
            zdvdt(jl,jk) = zdp(jl,jk) * &
                          (zmfuv(jl,ik)-zmfuv(jl,jk)+zmfdv(jl,ik)-zmfdv(jl,jk))
          end if
        end do
      else
        do jl = 1,klon
          if ( ldcum(jl) ) then
            zdudt(jl,jk) = -zdp(jl,jk)*(zmfuu(jl,jk)+zmfdu(jl,jk))
            zdvdt(jl,jk) = -zdp(jl,jk)*(zmfuv(jl,jk)+zmfdv(jl,jk))
          end if
        end do
      end if
    end do
!---------------------------------------------------------------------
!*  3.0        UPDATE TENDENCIES
!   -----------------
    do jk = ktopm2 , klev
      do jl = 1, klon
        if ( ldcum(jl) ) then
          ptenu(jl,jk) = ptenu(jl,jk) + zdudt(jl,jk)
          ptenv(jl,jk) = ptenv(jl,jk) + zdvdt(jl,jk)
        end if
      end do
    end do
!----------------------------------------------------------------------
    return
    end subroutine cududvn

!---------------------------------------------------------
!  level 3 souroutines
!--------------------------------------------------------
    subroutine cuctracer(klon,klev,ktrac,kctop,kdtop,           &
                       ldcum,lddraf,ztmst,paph,pmfu,pmfd,       &
                       pudrate,pddrate,pcen,ptenc)
    implicit none
    integer  klon,klev,ktrac
    integer  kctop(klon), kdtop(klon)
    logical  ldcum(klon), lddraf(klon)
    real(kind=kind_phys)   ztmst
    real(kind=kind_phys)   paph(klon,klev+1)
    real(kind=kind_phys)   pmfu(klon,klev)
    real(kind=kind_phys)   pmfd(klon,klev)
    real(kind=kind_phys)   pudrate(klon,klev)
    real(kind=kind_phys)   pddrate(klon,klev)
    real(kind=kind_phys)   pcen(klon,klev,ktrac)
    real(kind=kind_phys)   ptenc(klon,klev,ktrac)
    !----------------------------------------------------------------------
    integer ik , jk , jl , jn
    real(kind=kind_phys) zzp , zmfa , zerate , zposi
    ! ALLOCATABLE ARAYS
    real(kind=kind_phys) , dimension(:,:,:) , allocatable :: zcen , zcu , zcd , &
                               ztenc , zmfc
    real(kind=kind_phys) , dimension(:,:) , allocatable :: zdp
    logical , dimension(:,:) , allocatable :: llcumask , llcumbas
    !----------------------------------------------------------------------
    allocate (zcen(klon,klev,ktrac))  ! Half-level environmental values
    allocate (zcu(klon,klev,ktrac))   ! Updraft values
    allocate (zcd(klon,klev,ktrac))   ! Downdraft values
    allocate (ztenc(klon,klev,ktrac)) ! Tendency
    allocate (zmfc(klon,klev,ktrac))  ! Fluxes
    allocate (zdp(klon,klev))         ! Pressure difference
    allocate (llcumask(klon,klev))    ! Mask for convection
    ! Initialize Cumulus mask + some setups
    do jk = 2, klev
      do jl = 1, klon
        llcumask(jl,jk) = .false.
        if ( ldcum(jl) ) then
          zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk))
          if ( jk >= kctop(jl)-1 ) llcumask(jl,jk) = .true.
        end if
      end do
    end do
    !----------------------------------------------------------------------
    do jn = 1 , ktrac
      !*  1.0          DEFINE TRACERS AT HALF LEVELS
      !   -----------------------------
      do jk = 2 , klev
        ik = jk - 1
        do jl = 1, klon
          zcen(jl,jk,jn) = pcen(jl,jk,jn)
          zcd(jl,jk,jn) = pcen(jl,ik,jn)
          zcu(jl,jk,jn) = pcen(jl,ik,jn)
          zmfc(jl,jk,jn) = 0.
          ztenc(jl,jk,jn)= 0.
        end do
      end do

      do jl = 1, klon
        zcu(jl,klev,jn) = pcen(jl,klev,jn)
      end do
      !*  2.0          COMPUTE UPDRAFT VALUES
      !   ----------------------
      do jk = klev - 1 , 3 , -1
        ik = jk + 1
        do jl = 1, klon
          if ( llcumask(jl,jk) ) then
            zerate = pmfu(jl,jk) - pmfu(jl,ik) + pudrate(jl,jk)
            zmfa = 1./max(cmfcmin,pmfu(jl,jk))
            if ( jk >= kctop(jl) ) then
              zcu(jl,jk,jn) = (pmfu(jl,ik)*zcu(jl,ik,jn) + &
                zerate*pcen(jl,jk,jn)-pudrate(jl,jk)*zcu(jl,ik,jn))*zmfa
            end if
          end if
        end do
      end do
      !*  3.0          COMPUTE DOWNDRAFT VALUES
      !   ------------------------
      do jk = 3 , klev
        ik = jk - 1
        do jl = 1, klon
          if ( lddraf(jl) .and. jk == kdtop(jl) ) then
            ! Nota: in order to avoid final negative Tracer values at LFS
            ! the allowed value of ZCD depends on the jump in mass flux
            ! at the LFS
            zcd(jl,jk,jn) = 0.1*zcu(jl,jk,jn) + 0.9*pcen(jl,ik,jn)
          else if ( lddraf(jl).and.jk>kdtop(jl) ) then
            zerate = -pmfd(jl,jk) + pmfd(jl,ik) + pddrate(jl,jk)
            zmfa = 1./min(-cmfcmin,pmfd(jl,jk))
            zcd(jl,jk,jn) = (pmfd(jl,ik)*zcd(jl,ik,jn) - &
              zerate*pcen(jl,ik,jn)+pddrate(jl,jk)*zcd(jl,ik,jn))*zmfa
          end if
        end do
      end do
      ! In order to avoid negative Tracer at KLEV adjust ZCD
      jk = klev
      ik = jk - 1
      do jl = 1, klon
        if ( lddraf(jl) ) then
          zposi = -zdp(jl,jk) *(pmfu(jl,jk)*zcu(jl,jk,jn) + &
            pmfd(jl,jk)*zcd(jl,jk,jn)-(pmfu(jl,jk)+pmfd(jl,jk))*pcen(jl,ik,jn))
          if ( pcen(jl,jk,jn)+zposi*ztmst < 0. ) then
            zmfa = 1./min(-cmfcmin,pmfd(jl,jk))
            zcd(jl,jk,jn) = ((pmfu(jl,jk)+pmfd(jl,jk))*pcen(jl,ik,jn) - &
              pmfu(jl,jk)*zcu(jl,jk,jn)+pcen(jl,jk,jn) / &
              (ztmst*zdp(jl,jk)))*zmfa
          end if
        end if
      end do
    end do
    !----------------------------------------------------------------------
    do jn = 1 , ktrac
      !*  4.0          COMPUTE FLUXES
      !   --------------
      do jk = 2 , klev
        ik = jk - 1
        do jl = 1, klon
          if ( llcumask(jl,jk) ) then
            zmfa = pmfu(jl,jk) + pmfd(jl,jk)
            zmfc(jl,jk,jn) = pmfu(jl,jk)*zcu(jl,jk,jn) + &
              pmfd(jl,jk)*zcd(jl,jk,jn) - zmfa*zcen(jl,ik,jn)
          end if
        end do
      end do
      !*  5.0          COMPUTE TENDENCIES = RHS
      !   ------------------------
      do jk = 2 , klev - 1
        ik = jk + 1
        do jl = 1, klon
          if ( llcumask(jl,jk) ) then
            ztenc(jl,jk,jn) = zdp(jl,jk)*(zmfc(jl,ik,jn)-zmfc(jl,jk,jn))
          end if
        end do
      end do
      jk = klev
      do jl = 1, klon
        if ( ldcum(jl) ) ztenc(jl,jk,jn) = -zdp(jl,jk)*zmfc(jl,jk,jn)
      end do
    end do
      !*  6.0          UPDATE TENDENCIES
      !   -----------------
      do jn = 1, ktrac
        do jk = 2, klev
          do jl = 1, klon
            if ( llcumask(jl,jk) ) then
              ptenc(jl,jk,jn) = ptenc(jl,jk,jn)+ztenc(jl,jk,jn)
            end if
          end do
        end do
      end do
    !---------------------------------------------------------------------------
    deallocate (llcumask)
    deallocate (zdp)
    deallocate (zmfc)
    deallocate (ztenc)
    deallocate (zcd)
    deallocate (zcu)
    deallocate (zcen)
    end subroutine cuctracer

!---------------------------------------------------------
!  level 4 souroutines  
!--------------------------------------------------------
      subroutine cuadjtqn                 &                                           
     &    (klon, klev, kk, psp, pt, pq, ldflag,  kcall)                           
!          m.tiedtke         e.c.m.w.f.     12/89     
!          purpose.                                                                 
!          --------                                                                 
!          to produce t,q and l values for cloud ascent                             
                                                                                    
!          interface                                                                
!          ---------                                                                
!          this routine is called from subroutines:                                 
!              *cond*     (t and q at condensation level)                           
!              *cubase*   (t and q at condensation level)                           
!              *cuasc*    (t and q at cloud levels)                                 
!              *cuini*    (environmental t and qs values at half levels)            
!          input are unadjusted t and q values,                                     
!          it returns adjusted values of t and q                                    
                                                                                    
!     parameter     description                                   units             
!     ---------     -----------                                   -----             
!     input parameters (integer):                                                   
                                                                                    
!    *klon*         number of grid points per packet                                
!    *klev*         number of levels                                                
!    *kk*           level                                                           
!    *kcall*        defines calculation as                                          
!                      kcall=0  env. t and qs in*cuini*                             
!                      kcall=1  condensation in updrafts  (e.g. cubase, cuasc)      
!                      kcall=2  evaporation in downdrafts (e.g. cudlfs,cuddraf)     
!     input parameters (real(kind=kind_phys)):                                                      
                                                                                    
!    *psp*          pressure                                        pa              
                                                                                    
!     updated parameters (real(kind=kind_phys)):                                                    
                                                                                    
!    *pt*           temperature                                     k               
!    *pq*           specific humidity                             kg/kg             
!          externals                                                                
!          ---------                                                                
!          for condensation calculations.                                           
!          the tables are initialised in *suphec*.                                  
                                                                                    
!----------------------------------------------------------------------             
                                                                                    
      implicit none   
                 
      integer  klev,klon                                                                   
      real(kind=kind_phys)     pt(klon,klev),          pq(klon,klev),  &                             
     &         psp(klon)                                                            
      logical  ldflag(klon)                                                         
! local variables
      integer  jl,jk
      integer  isum,kcall,kk
      real(kind=kind_phys)     zqmax,zqsat,zcor,zqp,zcond,zcond1,zl,zi,zf
!----------------------------------------------------------------------             
!     1.           define constants                                                 
!                  ----------------                                                 
      zqmax=0.5  
                                                                                    
!     2.           calculate condensation and adjust t and q accordingly            
!                  -----------------------------------------------------            
                                                                                    
      if ( kcall == 1 ) then
      do jl = 1,klon
        if ( ldflag(jl) ) then
          zqp = 1./psp(jl)
          zl = 1./(pt(jl,kk)-c4les)
          zi = 1./(pt(jl,kk)-c4ies)
          zqsat = c2es*(foealfa(pt(jl,kk))*exp(c3les*(pt(jl,kk)-tmelt)*zl) + &
                (1.-foealfa(pt(jl,kk)))*exp(c3ies*(pt(jl,kk)-tmelt)*zi))
          zqsat = zqsat*zqp
          zqsat = min(0.5,zqsat)
          zcor = 1. - vtmpc1*zqsat
          zf = foealfa(pt(jl,kk))*r5alvcp*zl**2 + &
               (1.-foealfa(pt(jl,kk)))*r5alscp*zi**2
          zcond = (pq(jl,kk)*zcor**2-zqsat*zcor)/(zcor**2+zqsat*zf)
          if ( zcond > 0. ) then
            pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond
            pq(jl,kk) = pq(jl,kk) - zcond
            zl = 1./(pt(jl,kk)-c4les)
            zi = 1./(pt(jl,kk)-c4ies)
            zqsat = c2es*(foealfa(pt(jl,kk)) * &
              exp(c3les*(pt(jl,kk)-tmelt)*zl)+(1.-foealfa(pt(jl,kk))) * &
              exp(c3ies*(pt(jl,kk)-tmelt)*zi))
            zqsat = zqsat*zqp
            zqsat = min(0.5,zqsat)
            zcor = 1. - vtmpc1*zqsat
            zf = foealfa(pt(jl,kk))*r5alvcp*zl**2 + &
                 (1.-foealfa(pt(jl,kk)))*r5alscp*zi**2
            zcond1 = (pq(jl,kk)*zcor**2-zqsat*zcor)/(zcor**2+zqsat*zf)
            if ( abs(zcond) < 1.e-20 ) zcond1 = 0.
            pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
            pq(jl,kk) = pq(jl,kk) - zcond1
          end if
        end if
      end do
      elseif ( kcall == 2 ) then
      do jl = 1,klon
        if ( ldflag(jl) ) then
          zqp = 1./psp(jl)
          zqsat = foeewm(pt(jl,kk))*zqp
          zqsat = min(0.5,zqsat)
          zcor = 1./(1.-vtmpc1*zqsat)
          zqsat = zqsat*zcor
          zcond = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
          zcond = min(zcond,0.)
          pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond
          pq(jl,kk) = pq(jl,kk) - zcond
          zqsat = foeewm(pt(jl,kk))*zqp
          zqsat = min(0.5,zqsat)
          zcor = 1./(1.-vtmpc1*zqsat)
          zqsat = zqsat*zcor
          zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
          if ( abs(zcond) < 1.e-20 ) zcond1 = min(zcond1,0.)
          pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
          pq(jl,kk) = pq(jl,kk) - zcond1
        end if
      end do
      else if ( kcall == 0 ) then
      do jl = 1,klon
        zqp = 1./psp(jl)
        zqsat = foeewm(pt(jl,kk))*zqp
        zqsat = min(0.5,zqsat)
        zcor = 1./(1.-vtmpc1*zqsat)
        zqsat = zqsat*zcor
        zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
        pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
        pq(jl,kk) = pq(jl,kk) - zcond1
        zqsat = foeewm(pt(jl,kk))*zqp
        zqsat = min(0.5,zqsat)
        zcor = 1./(1.-vtmpc1*zqsat)
        zqsat = zqsat*zcor
        zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
        pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
        pq(jl,kk) = pq(jl,kk) - zcond1
      end do
      end if

      return
      end subroutine cuadjtqn
!---------------------------------------------------------
!  level 4 souroutines  
!--------------------------------------------------------
      subroutine cubasmcn &
     &    (klon,     klev,     klevm1,  kk,     pten,&
     &     pqen,     pqsen,    puen,    pven,   pverv,&
     &     pgeo,     pgeoh,    ldcum,   ktype,  klab,  plrain,&
     &     pmfu,     pmfub,    kcbot,   ptu,&
     &     pqu,      plu,      puu,     pvu,    pmfus,&
     &     pmfuq,    pmful,    pdmfup)
      implicit none
!      m.tiedtke         e.c.m.w.f.     12/89
!      c.zhang           iprc           05/2012
!***purpose.
!   --------
!          this routine calculates cloud base values
!          for midlevel convection
!***interface
!   ---------
!          this routine is called from *cuasc*.
!          input are environmental values t,q etc
!          it returns cloudbase values for midlevel convection
!***method.
!   -------
!          s. tiedtke (1989)
!***externals
!   ---------
!          none
! ----------------------------------------------------------------
      real(kind=kind_phys)     pten(klon,klev),        pqen(klon,klev),&
     &         puen(klon,klev),        pven(klon,klev),&
     &         pqsen(klon,klev),       pverv(klon,klev),&
     &         pgeo(klon,klev),        pgeoh(klon,klev+1)
      real(kind=kind_phys)     ptu(klon,klev),         pqu(klon,klev),&
     &         puu(klon,klev),         pvu(klon,klev),&
     &         plu(klon,klev),         pmfu(klon,klev),&
     &         pmfub(klon),   &         
     &         pmfus(klon,klev),       pmfuq(klon,klev),&
     &         pmful(klon,klev),       pdmfup(klon,klev),&
     &         plrain(klon,klev)
      integer  ktype(klon),            kcbot(klon),&
     &         klab(klon,klev)
      logical  ldcum(klon)
! local variabels
      integer  jl,kk,klev,klon,klevp1,klevm1
      real(kind=kind_phys)     zzzmb
!--------------------------------------------------------
!*    1.      calculate entrainment and detrainment rates
! -------------------------------------------------------
         do jl=1,klon
          if(.not.ldcum(jl) .and. klab(jl,kk+1).eq.0) then
            if(lmfmid .and. pqen(jl,kk) .gt. 0.80*pqsen(jl,kk).and. &
              pgeo(jl,kk)*zrg .gt. 5.0e2  .and.  &
     &        pgeo(jl,kk)*zrg .lt. 1.0e4 )  then
            ptu(jl,kk+1)=(cpd*pten(jl,kk)+pgeo(jl,kk)-pgeoh(jl,kk+1))&
     &                          *rcpd
            pqu(jl,kk+1)=pqen(jl,kk)
            plu(jl,kk+1)=0.
            zzzmb=max(cmfcmin,-pverv(jl,kk)*zrg)
            zzzmb=min(zzzmb,cmfcmax)
            pmfub(jl)=zzzmb
            pmfu(jl,kk+1)=pmfub(jl)
            pmfus(jl,kk+1)=pmfub(jl)*(cpd*ptu(jl,kk+1)+pgeoh(jl,kk+1))
            pmfuq(jl,kk+1)=pmfub(jl)*pqu(jl,kk+1)
            pmful(jl,kk+1)=0.
            pdmfup(jl,kk+1)=0.
            kcbot(jl)=kk
            klab(jl,kk+1)=1
            plrain(jl,kk+1)=0.0
            ktype(jl)=3
          end if
         end if
        end do
      return
      end subroutine cubasmcn
!
!---------------------------------------------------------
!  level 4 souroutines  
!---------------------------------------------------------
      subroutine cuentrn(klon,klev,kk,kcbot,ldcum,ldwork, &
                    pgeoh,pmfu,pdmfen,pdmfde)
       implicit none
       integer  klon,klev,kk
       integer  kcbot(klon)
       logical  ldcum(klon)
       logical  ldwork
       real(kind=kind_phys)  pgeoh(klon,klev+1)
       real(kind=kind_phys)  pmfu(klon,klev) 
       real(kind=kind_phys)  pdmfen(klon) 
       real(kind=kind_phys)  pdmfde(klon) 
       logical  llo1
       integer  jl
       real(kind=kind_phys)  zdz , zmf
       real(kind=kind_phys)  zentr(klon)
    !
    !* 1. CALCULATE ENTRAINMENT AND DETRAINMENT RATES
    ! -------------------------------------------
    if ( ldwork ) then
      do jl = 1,klon
        pdmfen(jl) = 0.
        pdmfde(jl) = 0.
        zentr(jl) = 0.
      end do
      !
      !*  1.1 SPECIFY ENTRAINMENT RATES
      !   -------------------------
      do jl = 1, klon
        if ( ldcum(jl) ) then
          zdz = (pgeoh(jl,kk)-pgeoh(jl,kk+1))*zrg
          zmf = pmfu(jl,kk+1)*zdz
          llo1 = kk < kcbot(jl)
          if ( llo1 ) then
            pdmfen(jl) = zentr(jl)*zmf
            pdmfde(jl) = 0.75e-4*zmf
          end if
        end if
      end do
    end if
    end subroutine cuentrn
!
!--------------------------------------------------------
! external functions
!------------------------------------------------------
      real(kind=kind_phys) function foealfa(tt)
!     foealfa is calculated to distinguish the three cases:
!
!                       foealfa=1            water phase
!                       foealfa=0            ice phase
!                       0 < foealfa < 1      mixed phase
!
!               input : tt = temperature
!
        implicit none
        real(kind=kind_phys) tt
         foealfa = min(1.,((max(rtice,min(rtwat,tt))-rtice) &
     &  /(rtwat-rtice))**2)

        return
      end function foealfa

      real(kind=kind_phys) function foelhm(tt)
        implicit none
        real(kind=kind_phys) tt
        foelhm = foealfa(tt)*alv + (1.-foealfa(tt))*als
      return
      end function foelhm

      real(kind=kind_phys) function foeewm(tt)
        implicit none
        real(kind=kind_phys) tt
        foeewm  = c2es * &
     &     (foealfa(tt)*exp(c3les*(tt-tmelt)/(tt-c4les))+ &
     &     (1.-foealfa(tt))*exp(c3ies*(tt-tmelt)/(tt-c4ies)))
      return
      end function foeewm

      real(kind=kind_phys) function foedem(tt)
        implicit none
        real(kind=kind_phys) tt
        foedem  = foealfa(tt)*r5alvcp*(1./(tt-c4les)**2)+ &
     &              (1.-foealfa(tt))*r5alscp*(1./(tt-c4ies)**2)
      return
      end function foedem

      real(kind=kind_phys) function foeldcpm(tt)
        implicit none
        real(kind=kind_phys) tt
        foeldcpm = foealfa(tt)*ralvdcp+ &
     &        (1.-foealfa(tt))*ralsdcp
      return
      end function  foeldcpm

      real(kind=kind_phys) function foeldcp(tt)
        implicit none
        real(kind=kind_phys) tt
        foeldcp = foedelta(tt)*ralvdcp + (1.-foedelta(tt))*ralsdcp
      end function foeldcp

      real(kind=kind_phys) function foedelta(tt)
        implicit none
        real(kind=kind_phys) tt
        foedelta = max(0.,sign(1.,tt-tmelt))
      end function foedelta

end module cu_ntiedtke
                                                                                         
